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ABSTRACT 



We use numerical simulations of critically-closed cold dark matter (CDM) models to 
study the effects of numerical resolution on observable quantities. We study simulations 
with up to 256^ particles using the particle-mesh (PM) method and with up to 144^ 
particles using the adaptive particle-particle-particle-mesh (P^M) method. Comparisons 
of galaxy halo distributions are made among the various simulations. We also compare 
distributions with observations and we explore methods for identifying halos, including 
a new algorithm that finds all particles within closed contours of the smoothed density 
field surrounding a peak. The simulated halos show more substructure than predicted 
by the Press-Schechter theory. We are able to rule out all O = 1 CDM models for linear 
amplitude erg ^ 0.5 because the simulations produce too many massive halos compared 
with the observations. The simulations also produce too many low mass halos. The 
distribution of halos characterized by their circular velocities for the P'^M simulations is 
in reasonable agreement with the observations for 150kms~^ ;S Vcirc ~ 350kms~^. 



Subject headings: cosmology: theory — dark matter — galaxies: clustering — galaxies: 
formation 



1. INTRODUCTION 



This paper is part of a two part series testing the cold dark matter (CDM) model of 
galaxy formation assuming a critically-closed universe, Q= 1. These papers focus on the 
formation and clustering of hales in cosmologically significant volumes of space (cubes 
of length ^ 50 Mpc on a side) with sufficient mass resolution and length resolution 
(force softening and box size) to resolve thousands of individual halos. The goal is not 
to study large scale structure ( ^ 200 Mpc; e.g. Park 1990). Rather, the goal is to 
study spatial and velocity statistics on scales ~ 1 — 10 Mpc using candidate galaxy halos 
identified in the nonlinear, evolved density field. A principal goal of both papers is to 
determine if there exists a linear normalization of the initial fiuctuation power spectrum 
(a free parameter in the theory) that satisfies observational constraints on galaxy masses, 
clustering, and velocities, and galaxy cluster multiplicity functions. 

The principle goal of this paper is to understand the properties of dark halos that 
form in cosmologically significant volumes of space in the CDM model. Specifically, we 
want to understand the sensitivity of halo formation and halo properties to numerical 
resolution. We identify which properties of halo formation {e.g. distributions of halo 
mass and circular velocity) are particularly sensitive to such parameters as box size, 
force resolution, mass resolution, and methods for identifying halos. 

Other workers have studied the formation of dark halos in the CDM scenario in 
volumes of space much greater than (100 Mpc)'^ by using approximate methods for iden- 
tifying galaxies as individual particles {e.g. Davis et al. 1985). Still others have studied 
volumes of space much smaller than (100 Mpc)^ with relatively high mass and force reso- 
lution {e.g. Frcnk et al. 1988). Small volumes of space do not contain long wavelengths in 
the initial conditions which may affect halo formation (studied in this paper) and which 
do affect clustering (Gelb & Bertschinger 1994, hereafter Paper II). The larger volumes 
of space simulated with relatively high numerical resolution presented in this paper also 
yield better statistics since more halos form than in smaller volumes. 

From our efforts, based on over one thousand IBM 3090 supercomputer- hours applied 
to more than a dozen large simulations, we gain insight into dynamic range by system- 
atically isolating various effects. We demonstrate which halo properties, if any, converge 
with increasing resolution up to practical limits using present-day supercomputers. These 
dynamic range studies are important for future workers who need to choose a particular 
set of simulation parameters for a particular problem in galaxy formation. 

By comparing the distribution of halo masses with estimates from observed galaxies, 
we show that the simulations produce too many massive halos. In Paper II, focusing on 
the spatial and velocity statistics of the halos, we consider the possibility that the overly 
massive halos represent clusters of galaxies (Katz & White 1993; Evrard, Summers, & 
Davis 1994). Because our simulations do not include gas dynamical dissipation, it is pos- 
sible that the dark matter halos we identify have clustering properties different from the 
luminous galaxies that would form if we properly simulated all of the physics of galaxy 
formation. To minimize the uncertainty caused by our lack of dissipative physics, we try 
to employ tests that should not depend strongly on the relation between dark halos and 



2 



luminous galaxies. For the same reason, in Paper II wc explore several different prescrip- 
tions for galaxy formation and we discuss cosmological N-body simulations employing 
gas dynamics {e.g. Cen & Ostriker 1992a,b; Katz, Hernquist, & Weinberg 1992). 

The N-body simulations follow the nonlinear gravitational clustering (in an expand- 
ing universe) of particles representing collisionless clouds of dark matter. The simulations 
utilize between 64^ (262144) and 256^ (16777216) particles in a universe with f2 = 1 and 
Hq = 50 kms"-*^ Mpc""*^. All distances are given in units Mpc rather than h~^Mpc. 
Most of the simulations are computed in cubes of length 51.2 Mpc on a side (box sizes 
are comoving). As we show in Paper II, this volume is too small to accurately measure 
galaxy clustering, although it allows one to resolve thousands of individual halos with 
hundreds to thousands of particles per typical Milky Way-sized halo. (We do, however, 
compute a few simulations in boxes of order 100 Mpc on a side in order to study galaxy 
clustering and small-scale peculiar velocities in Paper II.) 

Our simulations employ both the particle-mesh (PM) method (Hockney & Eastwood 
1982) and the adaptive particle-particle-particle-mesh (P^M) method (Couchman 1991). 
For a review of N-body methods in cosmology see Bertschinger (1991). Bertschinger 
& Gelb (1991) provide an overview of the numerical aspects of this work. Gelb (1992) 
provides many technical details and is the basis of these papers. 

In the remainder of this introduction we discuss briefly three key issues relevant for 
cosmological simulations of galaxy halo formation: force resolution, halo identification, 
and the normalization of the power spectrum. In § 2 we use the cumulative mass fraction 
of halos to study the effects of numerical resolution on halo formation, and we compare 
the simulations with the Press- Schechter (1974) theory. In § 3 we explore circular veloc- 
ity profiles and introduce observational data. In § 4 we compare the number of halos, 
characterized by their circular velocities, with observations. Separate sub-sections are 
included for high mass halos and for low mass halos. Conclusions and a summary are 
given in § 5. 

1.1. Force Resolution 

An important ingredient in N-body simulations is force resolution. We characterize 
the force softening in the simulations (with particle mass mpart) by the comoving pair 

separation r = R1/2 such that r"^ Fr / iGrn^.^^^) = 1/2, i.e. where the radial component 
of the force between two particles is half its Newtonian value. For the PM simulations 
Ri/2 ~ 1-4 grid cells (Gelb 1992 chapter 2). For P'^M simulations with a Plummer force 
law characterized by a softening e, i.e. with Fj. = Gm'^g^^^r / {r"^ + e^)'^/^, R1/2 ~ 1.305e. 
The shape of the PM softening is slightly different from a Plummer law, but in each case 
the appropriate force law (inverse square or Plummer) is matched accurately (to better 
than 2% rms) for r > 2f?i/2- There is additionally a small transverse component of the 
force due to PM grid anisotropics. Force errors are minimized using a suitable Green's 
function; see Bertschinger (1991), Gelb (1992), and Hockney & Eastwood (1982). 

For economy of notation and ease of reference we refer to the simulations as follows: 
CDMn(A'",L,i?i/2)- Following Gelb (1992), we number the simulations from n = 1 to 16. 
The numbers in parentheses indicate the following simulation parameters: 1) N particles. 
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2) a comoving box of length L Mpc on a side, and 3) a comoving force softening length of 
Ri/2 kpc. For example, CDM1(1283, 51. 2,280) uses 128^ particles, a (51.2 Mpc)^ box, and 
a comoving force softening length of 280 kpc. The two P'^M simulations discussed in these 
papers use R1/2 = 52 kpc comoving (e = 40 kpc) and R1/2 = 85 kpc comoving (e = 65 
kpc). The other simulations are low force resolution PM simulations with R1/2 > 190 
kpc comoving. 

We summarize the simulation parameters in Table 1. The entries are the following: 
simulation number, particle-mesh grid, particle mass, starting expansion factor, number 
of timesteps to (jg = 1, energy conservation relative to change in gravitational potential 
energy (see Gelb 1992 chapter 2), computer hours consumed, initial conditions identifier. 
Simulations with the same initial conditions identifier use equivalent sets of random 
numbers, i.e. they are generated from the same set of random numbers scaled to the 
appropriate power spectrum (see Gelb 1992 chapter 3). 

We use a time-centered leapfrog scheme (Hockney & Eastwood 1982) to advance the 
particles. All of the simulations are integrated using equal steps in expansion factor a, 
except CDM12, which uses equal steps in a°' with a = 0.5, as highlighted, for example, 
in the notable features column. 

All of the simulations use cloud-in-cell (CIC, see Hockney & Eastwood 1982) interpo- 
lation and a Holtzman (1989) CDM transfer function with 5% baryons, except CDM16, 
which uses triangular-shaped-cloud (TSC, see Hockney & Eastwood 1982) interpolation 
and a Bardeen et al. (1986, hereafter BBKS) transfer function. 

In order to avoid interference between the initial interparticle lattice and the particle- 
mesh grid (see Gelb 1992 chapter 2), we begin CDM6 with extra soft forces {i.e. we set 
the particle shape to be a linear sphere density profile with radius 77 = 5 grid cells, see 
Gelb 1992 Appendix I; then we set rj = 3.5 grid cells after the initial lattice disappears.) 
For CDM16, with 144^ particles, we use a 288^ grid (we use a 420^ grid after a = 0.7). 

High force resolution in a cosmologically significant box ( ^ 50 Mpc) is computation- 
ally challenging but can lead to significantly different results compared with low resolution 
simulations. One of our principle goals is to study the properties and clustering of re- 
solved halos, so we are forced to compromise mass and force resolution by using up to 
100 Mpc boxes. Other authors interested in the detailed properties of halos, and not 
clustering, have concentrated their efforts on very small box sizes. For example. Warren 
et al. (1991) used a tree code to simulate the formation of halos with very high particle 
number (1097921 particles) and very high force resolution (Plummer softening of 5 kpc 
proper) in a sphere of radius 5 Mpc. In another work, Dubinski and Carlberg (1991) 
studied CDM halo properties using a tree code with 32^ particles in a sphere of radius 
2.3 Mpc. The initial conditions were generated in a 8 Mpc box. The authors used an 
approximate treatment of tidal fields and a Plummer softening of 1.4 kpc. In the present 
paper the goal is to understand properties of halos evolved in larger boxes but with mass 
and force resolution significantly better than earlier efforts in boxes exceeding ~ 50 Mpc 
{e.g. Davis et al. 1985; White et al. 1987; Carlberg & Couchman 1989; Melott 1990; 
Park 1990). 

1.2. Halo Identification 
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The standard method for identifying halos from the evolved particle positions is to 
identify all particles within a given linking distance of each other (the friends-of-friends 
or FOF algorithm). We developed an alternative, novel procedure that identifies local 
density maxima in the smoothed, evolved density field: DENMAX (see Bertschinger 
& Gelb 1991; Gelb 1992 chapter 4). We first compute a static density field 5p/ p by 
interpolating the particles onto a grid. We then move the particles according to the 
equation 

f=v'-l, (1.1) 

dr p 

using a ficticious time variable r with 5p/ p held constant throughout the calculation. This 
equation describes a viscous fluid subject to a force proportional to the density gradient, 
in the limit of large damping. Every particle moves toward a density maximum where it 
comes to rest. All particles lying within closed density contour surfaces around a peak 
are pushed toward that peak. After the particles are sufficiently concentrated at density 
peaks, the particles are scooped up and their labels are recorded. A halo is composed of 
these particles with their original positions restored. The results of DENMAX depend on 
the degree of smoothing used to defined the density field 5p/ p. We use trilinear (cloud- 
in-cell) interpolation with a given grid {e.g., 512^ or smaller for sensitivity tests) to define 
the density field. 

After identifying halos, we remove the unbound particles, treating each halo in iso- 
lation. We compute the potential for each particle i, 0^, due to all Ny^ members of a 
halo: 

4'i = ^Hi^ij); = \ fj -fi\ . (1.2) 
3=1 

The potential is computed once and is fixed throughout the calculation. (For the P^M 
simulations we simply use the potential (/)(r) = —Gm/ (r^+e^)^/^. For the PM simulations 
we generate (/)(r) by integrating a force table generated by Monte Carlo sampling the PM 
force between pairs of particles.) We then iteratively remove unbound particles as follows. 
Wc compute the energy Ei = {l/2)m\vi — VcmP for each particle z, where v^ra is the 
mean velocity of the bound particles at any given stage. We then remove all particles 
with Ei > 0. The procedure is repeated, each time recomputing Vcm, until no more 
particles are removed. In all of the DENMAX analyses we remove the unbound particles. 
We have also identified halos using the FOF algorithm without the removal of unbound 
particles. 

1.3. Normalization of the Spectrum 

Most of the simulations are analyzed assuming three normalizations of the initial, 
linear CDM power spectrum (a free parameter in the theory). We define the normalization 
factor as using a tophat sphere of radius 8h~^ Mpc: 

oo 

al = J d^kPiin{k)W^ii{kR) ; R = Sh'^ Mpc (1.3) 
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with the tophat filter defined as 

WTnikR) = 



^ (sin kR — kR cos kR) 



(1.4) 



for comoving wavenumber k. The linear power spectrum of density fluctuations is 

Pii„(A;)= lim a-''P{k,ai) . (1.5) 

To define the CDM power spectrum we use the primordial scale-invariant spectrum 
modulated by the transfer function computed by BBKS or by Holtzman (1989) with 5% 
baryon fraction. The difference between the two is very small except at high wavenum- 
bers. We normalize the initial spectrum according to equation (1.3) with expansion factor 
a = 1 when erg = 1. We then scale the fluctuations to some early time a, using linear 
theory, i.e. P{k, ai) = a'fP{k, a = 1). 

We generally apply linear theory until the largest \Sp/p\ on the initial particle grid 
is imity. For the 144"^ particle simulation CDM16, however, linear theory is applied until 
the largest 3-dimensional displacement is 1 mean interparticle spacing, i.e. L/N^/^ for 
box size L and N particles. The Zel'dovich (1970) approximation is used to get particle 
positions and velocities at the end of the linear regime. The system is then evolved using 
the N-body code, with particle positions and velocities recorded at various expansion 
factors a = as- (By deflnition, erg oc a.) In most cases we study the models at erg = 0.5, 
0.7, and 1.0. In the literature, for example, ag = 0.4 is known as the h — 2.5 biased CDM 
model because of the assumption that galaxy density fluctuations are 2.5 times the mass 
density fluctuations. According to the linear biasing paradigm, h = l/crg. We do not 
adopt the linear biasing paradigm because we prefer to identify halos in the nonlinear, 
evolved mass distribution. Note that according to our prescription, the variance of halo 
numbers in d>h~^ Mpc spheres does not necessarily equal ag. The COBE measurement 
of microwave background anisotropy imply (for a scale-invariant spectrum of density 
perturbations and the standard CDM transfer function) erg ~ 1.1 (Wright et al. 1992; 
Efstathiou, Bond, & White 1992; Adams et al. 1993). 

2. DYNAMIC RANGE: CUMULATIVE MASS FRACTIONS 

In this section we discuss distributions of halos using the cumulative mass fraction 
(CMF). 

2.1. The CM¥ 

The CMF is deflned by: 

^ oo 

CMF(M) = ;^ J] ri{M') M', (2.1) 

M'=M 

where N is the total number of particles in the simulation, M is the mass (number 
of particles) of a halo, and n{M) is the number of halos containing M particles. By 
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definition, CMF(O) = 1, CMF(oo) = 0, and CMF(M) is a decreasing function of M. 
Note that the particle mass for N particles in a cube of comoving size L is 

m„„.(iV,L) =4.44x 10» Me {^) (5^)' • (2.2) 

The CMF gives the fraction of mass contained in halos more massive than M. Al- 
though the number and masses of large halos can fluctuate significantly from simulation 
to simulation, their contribution to the CMF gets averaged in the sum of equation (2.1). 
The smallest mass taken is typically 5 or more particles. The CMF has the advantage 
of summarizing in a nondimensional way all information about the mass function of ha- 
los. However, it has the disadvantage that halo masses are not easy to compare with 
observations. Also, because halos do not have sharp outer boundaries, the total mass 
of a given halo is often not a well-defined quantity. We address these problems later by 
applying a radial cut-off in order to compare with observations. Here the motives are 
purely theoretical in order to understand the effects of finite resolution. 

The first issue we study using the CMF is the difference between halos identified 
using DENMAX versus FOF. Fig. 1 shows the cumulative mass fraction versus mass 
for halos found in CDM1(1283,51.2,280) analyzed with DENMAX and FOF(/ = 0.1) 
and FOF(Z = 0.2), where I is the linking parameter in units of the mean interparticle 
spacing. The DENMAX masses include only the bound particles, while the FOF masses 
include all of the identified particles. The DENMAX CMFs lie between the FOF CMFs 
for / = 0.1 and / = 0.2. A smaller FOF linking parameter leads to smaller halos, but 
also to a smaller fraction of particles in halos. The reason for this is that FOF includes 
only particles such that the local over-density exceeds ~ 2l~^. DENMAX, however, 
gathers all particles around a peak, even those at lower density. FOF with I = 0.1 
dissolves low-density halos. If / is increased, then FOF merges halos together, increasing 
the maximum masses, even when the halos have distinct substructure (see Gelb 1992 
chapter 4, and Fig. 18 below). DENMAX avoids this problem: basically, any density 
concentration visible graphically will be found by DENMAX. (Indeed, graphical tests 
were first used to establish and test the algorithm.) Note that more than half of the 
particles are associated with some DENMAX halo, even at early times. This is the natural 
outcome of gravitational instability in a model with small-scale structure. Contrary to 
some expectations, most of the cold dark matter is not smoothly distributed. 

Although the differences in the CMF obtained using DENMAX and FOF are large, 
total halo masses are not measured in practice. It remains to be seen whether or not 
observable differences between DENMAX and FOF halos are large, and whether the 
results depend on the DENMAX grid or on I. DENMAX has a limitation stemming 
from the arbitrary choice of a density grid (512'^ for most of the analysis) or equivalently 
a smoothing scale for defining the density field. (Similarly, FOF has its own arbitrary 
parameter, /.) We explore these issues later. For now, our prejudice is to favor DENMAX 
because it does not suffer from the obvious defects of FOF, the dissolving of low-density 
halos and the merging of halos in high-density regions. We include FOF analysis only 
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for comparison with DENMAX because many authors use FOF {e.g. White et al. 1987; 
Carlberg, Couchman, & Thomas 1990; Brainerd & ViUumsen 1992). 

The lower panel in Fig. 1 shows the effect of the removal of unbound particles. The 
unbinding process systematically reduces the mass of the halos over the full range of 
masses, although the effects are largest for small masses. We find that the DENMAX 
results without the removal of unbound particles are in better agreement with I = 0.2 
FOF. However, the agreement is not exact; we show later that FOF occasionally links 
together visually distinct halos. Moreover, unbound particles are temporary members of 
the halos and therefore should not be included. 

Is there a significant simulation-to-simulation variation in the CMF? In Fig. 2 we 
show cumulative mass fractions for five simulations. They are all 128^ particle PM 
simulations (with R1/2 — 280 kpc comoving) computed in 51.2 Mpc boxes using different 
initial random numbers. There is very little scatter at the low mass end and there is 
considerable scatter at the high mass end. The fluctuations at the high mass end are 
due to small number statistics in these small volumes. We conclude that the CMF is not 
sensitive to simulation-to-simulation fluctuations except for rare massive halos. 

The next important issue is the effect of varying mass resolution and force resolution. 
In Fig. 3 we attempt to determine these effects by comparing four simulations in 51.2 
Mpc boxes which use initial conditions taken from an equivalent set of initial random 
numbers. (The same values are used for the initial Fourier transform of the density 
fluctuation fleld for all wavenumbers up to the Nyquist frequency for each cube. Thus, 
the initial conditions for = 128^ are identical to those for N = 64^ except that 
extra high-frequency power is present with the larger number of particles.) Mass and 
force resolution variations cause several effects that wc systematically separate out as we 
proceed. 

The N = 64^, R1/2 = 560 kpc comoving PM simulation fails to match up with 
the other simulations — this is not surprising considering that the force softening is so 
poor, larger than the size of many halos. The two very different simulations (the P"^M 
simulation with 64'^ particles and R1/2 — 52 kpc comoving versus the PM simulation with 
128^ particles and R1/2 — 280 kpc comoving) surprisingly yield very similar CMFs, but 
the harder forces in the P^M simulation actually give rise to halos with higher circular 
velocities, an important effect that is discussed in § 4. (We show as we proceed that the 
similarity of the CMF for these two simulations occurs because increased mass resolution 
and increased force resolution both increase the CMF.) The 256'^ particle simulation 
lies above the others due to the increase in mass resolution and the presence of more 
small-scale power in the initial conditions. 

2.2. DENMAX Resolution and Box Size 

We need to understand what happens if we vary the DENMAX grid when analyzing 
the same simulation. In Fig. 4 we show the results of several DENMAX analyses of the 
P^M simulation CDM12(643,51.2,52) at as = 0.5. We see that the DENMAX grid signif- 
icantly influences the CMF. This variation is analogous to the variation of the CMF with 
linking length I for the FOF algorithm (cf. Fig. 1). We demonstrate later, however, that 
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the circular velocities of the halos are less sensitive to the DENMAX grid — this is because 
circular velocities involve using a cut-off distance from the local density maximum. One 
effect arising from different DENMAX grids is the inclusion of distant particles into the 
halos. We demonstrate later that the DENMAX grid influences the break-up of massive 
halos when the grids are coarser than the force resolution of the simulation itself. Because 
of the density grid sensitivity of DENMAX, particularly for the total number of bound 
particles, we must compare the CMF from different simulations using the same effective 
DENMAX resolution. 

Are there significant differences in the CMF computed in boxes larger than 51.2 
Mpc? In Fig. 5 we show the CMF for two simulations computed in larger boxes (102.4 
Mpc and 100 Mpc). Since we also use a 512^ DENMAX grid for these simulations, the 
DENMAX resolution is only roughly half the resolution of the 51.2 Mpc simulations 
analyzed with a 512^ DENMAX grid. The difference is significant (cf. Kundic 1991). 
In order to separate out the effects due to larger waves in the initial conditions for 
the 100 Mpc boxes, compared with 51.2 Mpc boxes, we compare CDM16(144^,100,85) 
analyzed with a 512^ DENMAX grid with CDM12(643,51.2,52) analyzed with a 256^ 
DENMAX grid. (This is done at erg = 0.5 only.) The two simulations, CDM16 and 
CDM12, have comparable force resolution {R1/2 = 85 kpc comoving and 52 kpc comoving 
respectively) and comparable mass resolution (mpart = 2.3 x lO^'^ Mq and 3.5 x 10^*^ Mq 
respectively). The nearly perfect agreement between CDM16 (100 Mpc box) analyzed 
with 5123 DENMAX and CDM12 (51.2 Mpc box) analyzed with 256^ DENMAX, and the 
fact that the two simulations have comparable force and mass resolution, indicates that 
longer waves in the initial conditions do not significantly affect halo formation. (However, 
some of the longer waves have not gone nonlinear yet at cTg = 0.5.) This is encouraging 
because it means we can use the simulations in 51.2 Mpc boxes to understand halo 
properties. We will discover in Paper II, however, that the velocity dispersion of pairs of 
halos is significantly infiuenced by the different box sizes. 

To quantify the sensitivity of the CMF to DENMAX resolution, we measure the 
mass where the CMF equals 20%, denoted as M20. We choose 20% because larger val- 
ues are not well spanned by the various simulations and smaller values are more sensi- 
tive to the simulation-to-simulation variations of the massive halos. We compute, for 
CDM12(643, 51.2,52) at as = 0.5, the logarithmic slope A log M20/ A log D where D is the 
DENMAX grid spacing. In Fig. 4, comparing a 512^ DENMAX with a 256^ DENMAX 
analysis, we estimate AlogiQM2o ~ 12.65 — 13.33 = —0.68 and Alog^o-D = logio(l/2) 
so A log M20/ A log D ~ 2.27. Increasing D decreases the DENMAX resolution, thereby 
increasing the CMF. This is because a coarser DENMAX grid tends to pick out larger 
masses, i.e. it cannot resolve substructure. Comparing a 256^ DENMAX grid with a 
128^ DENMAX grid we find A logMso/A logi:> ^ 1.13. Comparing a 128^ DENMAX 
grid with a 64^ DENMAX grid we find A log M20/ A log L> ^ 0.57. We therefore see 
evidence for increasing amounts of substructure on smaller scales. Qualitatively similar 
behavior occurs with the FOF algorithm (cf. Fig. 1), where the linking parameter plays 
the role of the resolution scale. We demonstrate later that if we impose a radial cut on 
the DENMAX halos, as we do when we study circular velocities, the results are not as 
sensitive to resolution. 
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2.3. Small-scale Waves 

Fig. 6 is important for understanding the effect of varying the number of particles — 
particularly for separating out the fact that increasing the particle number not only 
increases the mass resolution, but it also probes smaller fluctuations in the initial power 
spectrum because of the higher Nyquist wavenumber cut-off. In a discrete system with N 
particles, the highest wavenumber represented, is {2n/ L){N^/^ /2) in each dimension. We 
show the results of a 512^ DENMAX analysis from three R1/2 = 280 kpc comoving PM 
simulations in 51.2 Mpc boxes which use equivalent initial conditions. The results are 
shown at three epochs for 128"^ and 64"^ particles. We also ran a simulation (CDM9) using 
128"^ particles, but the initial conditions are generated by interpolating the 64^ particle 
case to 128^ particles. Therefore, this simulation has the same mass resolution as the non- 
interpolated 128^ particle simulation (CDMl) but does not have the small-scale waves 
present in the non-interpolated simulation. 

Apart from the obvious increase in the CMF due to an increase in mass resolution 
(explored in greater detail below), we see in Fig. 6 the effect of the small-scale waves 
in the initial conditions — the non-interpolated 128"^ particle simulation has a higher 
value of the CMF at small mass relative to the interpolated 128^ particle case but not 
by much. Little, Weinberg, & Park (1991) studied the effect of the removal of high 
frequency waves in scale- free models. Using a PM simulation with 128"^ particles and 
P{k) oc k~^, they found that the nonlinear power spectrum in a simulation with initial 
power above kL/{2%) = 32 set to zero compared very well with the nonlinear power 
spectrum in a simulation with initial power above kL /{2tt) = 64 set to zero. Only small 
differences appeared on small scales, but further reductions in the initial cut-off frequency 
did produce large effects. 

2.4. Separation of Effects 

We now separate out the effects of mass and force resolution bearing in mind that 
1) we need to compare simulations in boxes of different sizes with the same effective 
DENMAX resolution; 2) the differences in the CMF arising from the inclusion of extra 
high and low frequency waves in the initial conditions are small; and 3) the simulation- 
to-simulation [i.e. different initial random numbers) differences in the CMF are small 
below about 10-^^ Mq. To separate out effects of resolution we re-examine Figs. 3 through 
6. 

We first demonstrate that higher mass resolution increases the CMF. If we examine 
Fig. 6 we see that the CMF is higher for the N = 128^ particle simulation than for the N = 
64^ particle simulation using the same force resolution (-R1/2 = 280 kpc comoving in a 51.2 
Mpc box) . Comparing 128"^ and 64'^ particle simulations, we find A log M20 / A log mpart ~ 
—0.56. The minus sign reflects the fact that if the particle mass increases, then the mass 
resolution oc 1/mpart decreases, and therefore M20 (or equivalently the CMF) decreases. 
The higher mass resolution simulations lead to a higher value of the CMF independent 
of force resolution. We also see this in Fig. 3 by comparing the 256*^ particle simulation 
(i?i/2 = 190 kpc comoving in a 51.2 Mpc box) with the 128*^ particle simulation {R1/2 — 
280 kpc comoving in a 51.2 Mpc box). The difference between R1/2 = 190 kpc comoving 
and Ri/2 = 280 comoving is shown later to have a nontrivial effect on the CMF. 
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For the 256^ particle simulation we find log^g M20 ~ 13.05 and for the 128^ particle 
simulation we find logj^g -^20 ~ 12.74. Therefore, we find AlogiVf2o/Alogmpart ~ 

-0.34. 

The effect on the CMF (logarithmic slope) is smaller going from 256^ to 128"^ particles 
(—0.34) compared with going from 128"^ to 64^ particles (—0.56), but it is not obvious if 
and when the results will converge. 

The fact that increased mass resolution continues to increase the CMF in the above 
comparisons warrants further investigation. Is this result still true when we impose a 
distance cut from the density peak? We re-analyze the three PM simulations (64^, 128^, 
and 256"^ particles) at ag = 0.5 imposing a distance cut of 300 kpc comoving from the 
density peak. The results are shown in the top panel of Fig. 7. In all three cases we 
do not remove the unbound particles from halos with raw masses (no cut in radius and 
no unbinding) exceeding 1.1 x IO^^Mq (location of vertical line; the transition mass) to 
be consistent with the analysis of the 256"^ particle simulation. (In all the analyses of 
the 256^ particle simulation CDM6 we do not remove the unbound particles from the 
massive halos, M > 1.1 x 10^'^ Mq, because it is computationally prohibitive.) However, 
the unbinding of the massive halos has a small effect on the CMF below the transition 
mass. To see this, we show in the middle panel of Fig. 7 at erg = 0.5 the CMF from 
the 64^ particle simulation and from the 128^ particle simulation with and without the 
unbinding of the massive halos. The effect is negligible just below the transition mass, 
and there is a slight increase in the CMF above the transition mass. 

By examining the top panel of Fig. 7 we find, for the 128"^ particle simulation versus 
the 64^ particle simulation, that A log M20 / A log TOpart ~ 

-0.23. This is less than -0.56, 
the result when we do not impose a cut of 300 kpc comoving from the density peak. 
The CMF itself changes considerably when we impose a distance cut from the density 
peak. However, we may adopt the position that particles at such great distances from the 
center of the halo should not be associated with estimated measurements of the mass of 
observed galaxy halos. The observed mass of individual galaxy halos at great distances, 
as opposed to dynamical properties inferred by the motions of stars and gas at small 
distances, is highly uncertain. By comparing the 256^ particle simulation with the 128^ 
particle simulation we find A log M20/ A log mpart ~ —0.22. Again this is less than —0.34, 
the result when we do not impose a cut of 300 kpc comoving from the density peak. 

Provided that we apply a cut in radius from the density peak, as we do when we 
characterize the halos by their circular velocities in the next sections, we see that the 
CMF is less sensitive to variations in mass resolution than when we do not impose a 
cut. We still do not see a convergence of the CMF with increasing mass resolution in 
Fig. 7ab. However, the 64^ particle simulation and the 128^ particle simulation both use 
Ri/2 = 280 kpc comoving. The 256^ particle simulation uses R1/2 = 190 kpc comoving. 
So next we correct for the difference in force resolution, but first we demonstrate that 
higher force resolution also increases the CMF. 

In order to see the effect of force resolution we compare simulations with similar mass 
resolution. In Fig. 5 we compare CDM12(64^,51.2,52) analyzed with a 256^ DENMAX 
grid and CDM11(1283, 102.4,560) analyzed with a 512^ DENMAX grid. The mass reso- 
lution and the DENMAX grid resolution are equivalent since the CDMll box has eight 
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times the volume of the CDM12 box. We see that the higher force resolution simulation 
yields a higher value of the CMF. We find that M20 ~ 13.33 for the high resolution simu- 
lation and 13.08 for the low resolution simulation. If we characterize the force resolution 
by -Ri /2 , then we find that A log M20 / A log Ri /2 ~ —0.24. This number should be treated 
with caution since we are comparing simulations with R1/2 = 52 kpc comoving versus 
Ri/2 — 560 kpc comoving — this is a wide range and DENMAX behaves unreliably in very 
low resolution simulations. We do not have two P"^M CDM simulations with comparable 
mass resolution but with significantly different Plummer softenings. 

The increase in the CMF for higher force resolution simulations is verified by com- 
paring CDM12(643,51.2,52) with CDM8(643, 51.2,560) in Fig. 3, but again the force res- 
olution in CDM8 is extremely poor. 

As a final comparison of force resolution effects, we compare CDM12(64^,51.2,52) 
with CDM7(64^51.2,280) both analyzed with a 512^ DENMAX grid. We find 
A log M20/A log i?i/2 ~ —0.56. The range of force softenings in this comparison is still 
large but at least R1/2 — 280 kpc comoving is more reasonable than 560 kpc comoving. 
In subsequent sections we compare the halos characterized by their circular velocities 
and particular attention is paid to force resolution comparing results for PM versus P^M 
simulations. So we return to force resolution then. 

As a final test of the convergence of the CMF with increasing mass resolution, we first 
use the above force resolution analysis to estimate the effect on the CMF from a R1/2 = 
280 kpc comoving PM simulation versus a R1/2 = 190 kpc comoving PM simulation. 
To do so, we compare the 64^ particle, R1/2 = 52 kpc comoving P^M simulation with 
the 64^ particle, R1/2 = 280 kpc comoving PM simulation (both computed in 51.2 Mpc 
boxes and analyzed with a 512^ DENMAX grid) imposing a 300 kpc comoving cut from 
the density peak of the halos. The logarithmic slope is A logM2o/A logi?i/2 ~ 0.74. If 
we multiply —0.74 by Alogio-Ri/2 — —0.18, i.e. the difference between the softening 
of the Ri/2 = 190 kpc comoving PM simulation and the R1/2 = 280 kpc comoving PM 
simulation, we get Alog^Q M20 ~ 0.13. Therefore, we can estimate that the 256^ particle 
Ri/2 = 190 kpc comoving PM simulation (with a 300 kpc comoving distance cut) would 
have log]^Q M20 ~ 11.91 — 0.13 = 11.78 if it were computed using a R1/2 = 280 kpc 
comoving PM simulation. 

Now if we compare the re-scaled (to R1/2 = 280 kpc comoving) 256"^ particle result 
with the Ri/2 = 280 kpc comoving 128^ particle PM simulation, all with a 300 kpc co- 
moving distance cut, we get A log M20/ A log mpart ~ —0.11 compared with the old value 
of —0.22. This is encouraging because this logarithmic slope, —0.11, is still better than 
the logarithmic slope —0.23 computed earlier by comparing a 128^ particle simulation 
with a 64^ particle simulation. Graphically (as depicted in Fig. 7c) this corresponds to 
moving the CMF for the 256^ particle simulation in the top panel of Fig. 6 0.13 units to 
the left. 

We now see that the agreement between the 128^ particle case and the 256^ particle 
case is much better (Fig. 7c). There is still a slight increase in the CMF on small mass 
scales. However, this is consistent with the fact that the 256"^ particle simulation has more 
small scale power in the initial conditions compared with the 128^ particles simulation. 
This effect was demonstrated earlier. 
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Using simulations analyzed with the same effective DENMAX resolution, we found 
the following: 1) Higher mass resolution leads to larger values of the CMF independent 
of force resolution. The effect is smaller when we impose a distance cut from the density 
peak of the halos. The difference between the 4.4 x lO^M© and the 3.5 x 10^^ Mq simu- 
lation (using a distance cut of 300 kpc comoving from the density peaks of the halos) is 
small, A log Af2o/ A log rnpart ~ —0.23. The difference between the 5.5 x lO^M© simula- 
tion and the 4.4 x IO^Mq simulation (using a 300 kpc comoving cut and correcting for 
the difference in force softening) is A log M20 / A log mpart ~ —0.11. The difference has 
decreased in the very high mass resolution simulation indicating that convergence of the 
CMF with mass resolution is plausible. 2) Higher force resolution leads to larger values 
of the CMF independent of mass resolution. We examine the effects on the formation 
of halos arising from different force resolution P^M simulations in subsequent sections. 
3) Longer waves in the initial conditions (100 Mpc box versus a 51.2 Mpc box) do not 
significantly affect the CMF. 4) Smaller waves in the initial conditions (64^ particle initial 
conditions interpolated to 128^ particles versus true 128^ particle initial conditions) do 
not significantly affect the CMF, aside from a small effect on small mass scales. 5) Larger 
DENMAX grids better resolve substructure; this lowers the CMF. The results are sensi- 
tive to the different DENMAX grids so it is important to compare CMFs using the same 
effective DENMAX resolution. However, we show later that the results are less sensitive 
when we compute circular velocities which are what we use to compare the simulated 
halos with the observations. 



2.5. Press-Schechter Theory 

As a final application of the CMF, we compare the simulations with the predictions 
of the Press-Schechter theory (Press & Schechter 1974). The Press-Schechter formalism 
estimates the fraction of mass in bound halos with masses > M to be the fraction of the 
mass whose linear density, averaged over a scale M, exceeds 5^. 



P{M) = erfc 



2V2ao(M) 



(2.3) 



where erfc is the complementary error function. One may regard 5c as a free parameter, 
although it is often taken to equal the critical over-density for uniform spherical collapse 
in an Einstein-de Sitter universe, 5c — 1.68. The rms density (Jo(M) is computed from 
the linear power spectrum, smoothed with an appropriate filter (window function). We 
use either a Gaussian window function, W{k.,Rf) = exp(— 0.5a;^), or a tophat window 
function, W{k.,Rf) = 3(sinx — xcosx)/x^, where in both kRf. The generalized 

spectral moments (to be used below) are defined as follows: 



00 

al{M) = J A7rk^P{k) W^{k, Rf)P''dk . (2.4) 


For a Gaussian window function, the smoothing radius Rf is related to the mass as 
follows: M = (27r)^/^po-R/- For a tophat window function, M = (47r/3)i?^. 
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Press &; Schechter (1974) estimated the mass function of bound halos as 
n{M)dlnM = 2po{dP/dM)d\ia.M , where po is the comoving background mass density. 
The factor of two is needed for normahzation, but has since been derived analytically by 
Bond et al. (1991). The final result is: 



n(M)dlnM=(-V^'^^^ 



d In o"o 



dlnM 



exp 



Sr. 



2\ao{M) 



dlnM . 



(2.5) 



We convert equation (2.5) into a CMF using 

oo 

CMF(M) = — I n{M')M'd\nM' . 
Po J 



(2.6) 



InM 



We evaluate equation (2.6) using erg = 0.5, 0.7, and 1.0 linear normalizations of 
the BBKS CDM power spectrum. We try 5c = 1.44 {e.g. Carlberg & Couchman 1989), 
Sc = 1.68 {e.g. Efstathiou et al. 1988; Brainerd & Villumsen 1992), and Sc = 2.0, for both 
a gaussian and a tophat window function. Theoretical predictions of the Press- Schechter 
theory are compared with CMFs measured from the high resolution N-body simulations 
CDM12(643,51.2,52; particle mass 3.5 x lO^^M©) and CDM16(1443, 100,85; particle mass 
2.3 X IO^^Mq) in Fig. 8. 

First we consider the halos identified according to the FOF algorithm with a linking 
parameter / = 0.2. Fig. 8a shows that the two simulations, at three different epochs, yield 
reasonably good agreement with the Press-Schechter predictions for a tophat window 
function with Sc — 2.0. Only slightly worse agreement obtains with a gaussian window 
function with Sc = 1.68. Note that the simulated mass distributions are broader than 
predicted. The high mass tails of the distributions actually match very well the Press- 
Schechter predictions for a tophat window function with Sc = 1.68, but there are fewer 
low mass halos than predicted. Evidently this is because they are subsumed into more 
massive halos, at least with the FOF recipe, with greater efficiency than implied by 
the analytical model. Although the agreement with the Press-Schechter theory is not 
perfect, the errors do not grow with epoch; the analytical theory appears to give the 
correct scaling of masses as the clustering strength increases. Our result here differs from 
that of Brainerd & Villumsen (1992), who found the departures growing as clustering 
progresses. 

Fig. 8b shows similar results for a FOF linking length I = 0.1. Now Sc must be 
increased (from 1.68 to 2.0 for the gaussian window function) to account for the smaller 
masses of the halos defined at a higher over-density. However, the agreement at small 
masses is significantly worsened. 

Fig. 8c shows CMFs for CDM12(643,51.2,52) computed using DENMAX, compared 
with Press-Schechter theory for a gaussian window function with Sc — 2.0. The top set of 
data points (filled circles and solid curves) are for raw DENMAX masses, with no removal 
of unbound particles (which would decrease the CMF by about 10%) and with no radial 
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cut. The bottom set (dashed curves and crosses) have excluded unbound particles and 
those beyond a comoving radius of 200 kpc from the peak. There are several important 
things to notice. First, at early epochs, the raw DENMAX CMF agrees well with the 
Press-Schechter theory. At high masses the DENMAX distributions are similar to those 
obtained using FOF with I = 0.1 while at low masses they match the I = 0.2 case better. 
DENMAX breaks up the more massive clumps found with I = 0.2 while preserving the 
subclumps as individual halos. 

Second, as clustering increases, the CMF grows less rapidly than the Press-Schcchter 
prediction. This effect appears to be due to the ability of DENMAX to find substructure 
in halos merged by FOF. Thus, although we disagree with Brainerd & Villumsen (1992) 
about the results from FOF, we agree that the actual halo mass distribution grows less 
rapidly than predicted by Press-Schechter theory. The agreement could be improved if 
5c were to grow with epoch. In fact, at very early epochs (when there are fewer than 
100 particles per group) the fit to the simulations is good with a gaussian window and 
6c = 1.68. 

The third point to note from Fig. 8c is that the radial truncation of the halos makes 
a big difference in the masses. Thus, the halos are very extended, a point that we will 
demonstrate more clearly later. 

In summary, halo mass functions depend on how the halos are defined. Earlier 
workers {e.g., Efstathiou et al. 1988; Carlberg & Couchman 1989) found good agreement 
between the Press-Schechter theory and simulations. However, the simulations were ana- 
lyzed with a low resolution group finder, FOF (^=0.2), and the halos contained relatively 
few particles. Our results agree with this work, but show further that the Press-Schechter 
theory does not match well the CMF when higher resolution is used to identify halos made 
of thousands of particles. The disagreement is in the sense that the simulated halos are 
less massive than predicted. This occurs not because large halos have failed to collapse. 
Rather, merging does not immediately erase the substructure in large halos, contrary to 
the assumptions made in the Press-Schechter theory. 

3. DISTRIBUTIONS OF HALOS: BACKGROUND 

3.1. The Schechter Luminosity Function 

We need to define physically motivated catalogs of halos in order to understand 
further the effects of dynamic range on halo formation and in order to compare the 
simulations with the observations. Total bound mass, as in the previous section, is only 
one way to characterize the halos. We can also ask how much mass is contained within 
a specified radius. This is equivalent to specifying Vdrc = {GM/Ry/'^. Empirically, Vdrc 
is found to be nearly independent of R and to correlate well with optical luminosity. We 
will use these correlations — the TuUy-Fisher (1977) relationship for spiral galaxies and 
the Faber-Jackson (1976) relationship for elliptical galaxies — to assign a luminosity to 
each halo. 

Observations of spiral galaxies are measured in terms of their circular velocity and 
observations of elliptical galaxies are measured in terms of their average central radial 
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velocity dispersion. (Technically, the elliptical observations are luminosity weighted mea- 
surements of radial velocities along the line-of-sight.) 

We realize that we cannot adequately relate internal velocity dispersions of dark 
matter to velocity dispersions of centrally concentrated stars. Nevertheless, we define 
the quantities cJi and 0"^- (cr^ is closer to what the observers measure) from the simulated 
galaxy halos as follows: 




where represents the number of bound particles within a distance R from the local 
density maximum and fi is the unit vector from the local density maximum to particle 
i. We do not attempt to distinguish the simulated halos as spirals or ellipticals; rather, 
we characterize all of the simulated halos in terms of their circular velocities. 

Because the velocity dispersion tensor is radially anisotropic we find that ai is typi- 
cally ~ 20% lower than ar- We study both quantities, using various cut-off radii (typically 
a few hundred kpc comoving), when comparing the velocity dispersions of massive simu- 
lated halos (perhaps associated with elliptical galaxies) with observations of the velocity 
dispersions of centrally concentrated stars. In order to test if either cxi or Cr is a useful 
statistic, and because the stars are in orbits with smaller apapses than the dark matter, 
we use a crude, linear scaling law (derived from observations of M87) as discussed in 
greater detail in § 4.4. 

For = 1 and h = 1/2, the circular velocity, for an assumed spherical halo, as a 
function of total particle number in the simulation, N, and the comoving box size in Mpc, 
L, is: 

..„(.) = ....s-y^(l|!)"^(^)''^(^)''^ , (3.) 

where Nc{R) is the total number of bound particles within a comoving distance R from 
the smoothed density maximum found by DENMAX. The present epoch is a = oq = 
erg. In most of the figures, we assume that a = ao and we consider different possible 
normalizations by varying cxg — do- In one case below (Fig. 16), we fix ao and look at the 
evolution of halos for different a. In all cases, we take to be a comoving radius (i.e., 
a proper radius at a = ao) and we use the proper radius Ra/ao in the denominator. To 
get circular velocities measured at a fixed comoving radius, we set a = ao. 

In order to compute the observed distribution of galaxies as a function of Kirc, i-^- 
-^(Kirc)AV'circ, wc assumc a Schechter (1976) luminosity function with the form 

^C)dC = $* exp {-C/C,){C/C.)'^d{C/C,) , (3.3) 

where <^{C)dC is the density of galaxies in the luminosity range Cio C + dC. We convert 
equation (3.3) into counts of halos in a (51.2 Mpc)^ comoving volume as a standard 
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reference for all of the simulations in bins of V^irc using a relationship for L = C{Vcirc)- 
We also use blue magnitudes and selected values of and (both assuming h = 1/2), 
and a value of a. 

We define the distribution of halos, or number of halos binned by V^iro as 

iV(Feirc)AFeirc = (^^il^i^E^^ N {V,,rc) AV,i,, , (3.4) 

where A'^(Vcirc) is the number of halos found in the simulation with circular velocities in 
the range 

Kirc ± AV;irc/2 with AFcirc = 25kms-^ The factor (51.2 Mpc/L)^ is used to 
scale all of the results to comoving volumes (51.2 Mpc)^ for comparison. 

We compute the corresponding mean number of galaxies from the observations as fol- 
lows, assuming C is related to Mbt and Mb^ — /(Kirc) for some function / given below: 

iVschechter(Kirc)AKirc = (51.2 Mpcf J $(x) dx 

=/(Fcirc + AFeirc/2) 
Mg^ = /(Kirc - AKirc/2) 

^{x[Mb^]) = $*x"exp(-a;) ; x = iQi^B^-^B^)/'^-^ ^ ^ . (3.5) 

We use the central values of parameters found by Efstathiou, Ellis, & Peterson (1988): 
= (1.56 ± 0.34) X 10-2/^3 Mpc"^ M*^ = -19.68 ± 0.10 - 2.51oglo/^"^ and a = 
-1.07 ±0.05. 

For the function /(V^circ) for spiral galaxies we use the blue TuUy- Fisher relation from 
Pierce & TuUy (1988): 

/spiral (K:irc) = Mb^ = -6.86 log^Q (2Fcirc) ' 2.27 + 5 logiQ (50/85) + 0.569 . (3.6) 

The term 5 logiQ(50/85) is used to convert from a Hubble constant of 85 kms~^ Mpc~^ 
to 50 kms~^ Mpc~^. The term 0.569 is used to correct for random inclinations following 
TuUy & Fouque (1985). 

For the function /(V^circ) for elliptical galaxies we use the Faber-Jackson relation 
from our fit (unpublished) to elliptical data of Faber et al. (1989), assuming a Hubble 
constant of 50 kms~^ Mpc~^: 

/elliptical (Kirc) = Mb^ = -6.6364 logio((Ti) - 5.884 , (3.7) 
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where we relate ui to V^irc using: 



= . (3.8) 



The factor F, discussed in the next section, is estimated from the simulations. This use 
of o"!, however, is an oversimplification (mostly affecting high Vdrc) for reasons discussed 
earlier. Again, we re-examine the high mass halos in detail in § 4.4, where we use ui, ar, 
and a linear scaling law derived from M87. 

The final ingredient is to assume that 70% of the galaxies are spirals and 30% are 
ellipticals. This is also the assumption used by Frenk et al. (1988). In other words, we 
add together the results for spirals using equation (3.6) to relate circular velocities to 
absolute magnitudes and weighting equation (3.5) by 0.7 with the results for ellipticals 
using equation (3.7) to relate circular velocities to absolute magnitudes and weighting 
equation (3.5) by 0.3. Dressier (1980), however, found a higher concentration of ellipticals 
in rich clusters compared with lower density regions. Postman & Geller (1984) found for 
the CfA survey that 1) the relative numbers of galaxies are 65% spirals, 23% SO's, and 
12% ellipticals and 2) there is a dramatic increase in the relative number of spirals in the 
field compared with dense regions. These percentages can alter the estimates at the high 
mass end. 



3.2. (7i versus VcWc 

The factor F in equation (3.8) is measured empirically from the simulation 
CDM16(144^100,85) using DENMAX halos analyzed with a 512^ grid. White et al. 

(1987) used F = 1 (in our notation) but the same authors used F = 1.1 in Frenk et al. 

(1988) . 

Because the Plummer softening in CDM16 is e = 65 kpc comoving (or R1/2 = 85 kpc) 
we cannot directly determine o'i{R) or V"circ(-R) at the distances where optical observations 
of real galaxies are made. Optical observations of central velocity dispersions of large 
elliptical galaxies are made on scales of a few kpc to ~ 6 kpc (see Franx, lUingworth, 
& Heckman 1989). Optical observations of circular velocities of large spiral galaxies are 
made out to ~ 10 kpc. Rubin et al. (1985) studied 16 large spiral galaxies where they 
could measure velocities out to large radii. The average maximum distance for which 
they made measurements was 16.4 kpc and the maximum distance for the 16 galaxies 
was 51.2 kpc. We consider the limitations arising from our measurements at large radii 
as we proceed. 

What is a value of R for computing ai where the results are independent of R? In 
Fig. 9 (top panel) we show cxi evaluated using R ~ 100 kpc comoving versus R = 200 
kpc comoving. The slight increase in ai for R = 200 kpc comoving versus R = 100 kpc 
comoving (top panel) indicates that contributions from particles at large separations are 
still important for the most massive halos. We find this trend to be larger when comparing 
results from i? = 50 kpc comoving versus i? = 100 kpc comoving, indicating that i? = 100 
kpc comoving is too small. We find this trend to be small when comparing results from 
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R = 200 kpc comoving versus R = 300 kpc comoving, indicating that R = 200 kpc 
comoving is adequate. We find similar results at ag = 0.7 and erg = 1.0. 

What is a value of R for computing V^irc where the results are independent of R7 
In Fig. 9 we also show computations for V"circ using R = 100 kpc comoving, 200 kpc 
comoving, and 300 kpc comoving. The results indicate that R = 200 kpc comoving is 
acceptable (bottom panel) while R = 100 kpc comoving again is too small (middle panel). 

What is an empirical value of F in equation (3.8)? In Fig. 10 we show V"cii-c/c"i versus 
V^irc (all computed with R — 200 kpc comoving) at ag — 0.5, 0.7, and a=1.0 for halos 
from CDM16(1443, 100,85). The solid lines indicate F =1 {i.e. Fcirc/(^i = VS) and the 
dotted lines indicate F = 1.1 (i.e. Vdrc/cn = v^/1.1). There is less scatter for high 
values of Vcirc versus low values of Vdrc- The factor F affects the conversion of Vdrc to 
cTi for ellipticals. Ellipticals dominate at the high mass end where, at erg = 0.5, F = 1.1 
works slightly better than F — 1. However, when we show the observed number of halos 
in § 4 we use both F = 1 and F = 1.1 — the latter yields fewer bright halos since it 
effectively raises ai for a given V^irc implying a brighter elliptical galaxy (or equivalently, 
F = 1.1 effectively assigns a smaller stellar velocity dispersion for a given ui). 

In summary, we compute the number of halos scaled to (51.2 Mpc)^ comoving 
volumes assuming a Schechter luminosity function ($* — 1.56 x 10~^h^ Mpc~"^ and 
M*^ = -19.68 - 2.51ogio h'^ with h = 1/2 and a = -1.07 ). We assume 70% of the 
galaxies are spirals with a TuUy-Fisher relation given by equation (3.6) and 30% of the 
galaxies are ellipticals with a Faber- Jackson relation given by equation (3.7). We convert 
elliptical measurements in terms of ai to Vch-c using equation (3.8) for both F = 1 and 
F = 1.1. For the most massive halos F = 1 is adequate except at us = 0.5 where F = 1.1 
is slightly better. 
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3.3. Circular Velocity Profiles 

We now examine circular velocity profiles and 1-dimensional velocity dispersion pro- 
files from CDM12(64^,51.2,52) in Fig. 11. We extract several facts from these plots. First, 
the circular velocities are very flat for nearly all the halos, except for a few massive ones, 
for R ^ 150 kpc comoving. We also see that the circular velocities are flat for many of 
the midsize halos down to about 80 kpc comoving (twice the Plummer softening length) . 
We conclude that 150 kpc comoving is a good place to characterize the circular velocities 
for this simulation at all epochs. We arc making a signiflcant error only for a handful of 
the most massive halos. The rising circular velocities for the most massive halos are not 
an artifact of softening (see § 4.4). 

We also study CDMl (i?i/2= 280 kpc comoving) and CDM16 (i?i/2 = 85 kpc co- 
moving) and we find that R = 300 kpc comoving is suitable for CDMl and R = 200 
kpc comoving is suitable for CDM16. In most cases, we use R = 150 kpc comoving 
for the Ri/2 = 52 kpc comoving simulation; R = 200 kpc comoving for the R1/2 = 85 
kpc comoving simulation; and R = 300 kpc comoving for the R1/2 = 280 kpc comoving 
simulation. These values are chosen where Vcirc is flat for nearly all of the halos. 

We compare our circular velocity proflles to simulations by others with much higher 
force resolution. The fl = 1 CDM simulations of Dubinski & Carlberg (1991) used 33000 
particles in a 2.3 Mpc radius sphere. The typical particle mass is 1.2 x 1O^M0 and the 
Plummer softening is 1.4 kpc. In their figure 4 they show several circular velocity profiles 
with halos that have maximum circular velocities of about 290 kms~^. These halos have 
flat circular velocities between about 10 kpc and 60 kpc (there is a very slight decrease 
over this range) . The circular velocities rise on a scale a few times the Plummer softening 
length as we also find in our simulations. 

The simulations of Warren et al. (1991) used roughly a million particles in a 5 Mpc 
radius sphere and a Plummer softening of 5 kpc. For circular velocity profiles that have 
maximum circular velocities of about 150kms~^ they foimd that the profiles are rising 
out to a distance of about 30 kpc — again, several Plummer lengths. They also found 
that their profiles are falling typically beyond a distance of about 40 kpc. In our P^M 
simulations we do not find falling circular velocity profiles until a distance of about 100 
kpc. One reason for this disrepancy is that Warren et al. did not not use a CDM power 
spectrum — rather, they used P{k) oc /c on large scales with a sharp transition, at 1.5 
Mpc, to P{k) oc on small scales. The behavior of circular velocities is a function of 
the effective index in the initial power spectrum (Hoffman &; Shaham 1985). 

We assume that if we had used force softening below the typical ~10 kpc observed 
scale, as in the simulations of Dubinski & Carlberg (1991) and Warren et al. (1991), that 
our circular velocity profiles might remain fiat down to these scales. Therefore, we do not 
expect to make a significant error by estimating V^circ using R ^ 150 kpc comoving. We 
cannot use circular velocities to characterize the most massive halos because observational 
data for massive halos are based on velocity dispersions, not circular velocities. We 
explore the properties of cri below. 

There is cause for concern when using simulations with force softening far beyond a 
few kpc — the scale beyond which most spiral galaxies are observed to have fiat rotation 
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curves. Arc wc able to adequately resolve individual halos? We comment on several 
issues related to this question. 1) Using ~ (1/$*)^/^, and $* = 1.95 x 10""^ Mpc~^, 
we find the mean spacing between bright galaxies is ~ 8 Mpc. This is much greater 
than our fiducial radius ~ 200 kpc. Observations show flat rotation curves "as far as the 
eye can see" for most spiral galaxies (Rubin et al. 1985). It is therefore possible that 
real galaxies have flat rotation curves beyond 200 kpc. 2) The mean galaxy spacing is 
much smaller in rich clusters. It is possible that some of our massive halos are mergers 
where dissipative effects might allow many galaxies to survive in a single halo (White & 
Rees 1978; Katz & White 1993). In Paper II we break up these systems using various 
methods in an attempt to estimate the effects on clustering. In this paper, however, we 
consider the massive halos at face value and we examine the implications for CDM in § 
4. 3) We compare results from the e = 40 kpc comoving simulation with results from the 
e = 65 kpc comoving simulation and we ask if the distributions of halos are significantly 
different. 

Last, we show (Ti{R) for the e = 40 kpc comoving simulation CDM12(64^,51.2,52) 
in Fig. 12. (These are the same halos shown in Fig. 11.) The first thing we notice is 
that the profiles are very fiat down to about 40 kpc comoving, the Plummer softening 
scale. (On smaller scales we are limited by both force and mass resolution.) Also, the 
profiles for the most massive halos are fiat down to typically 100 kpc comoving. In § 4.4 
we use ar computed at large radii (similar to ui) to compare simulated halos with the 
observations. 

4. DISTRIBUTIONS OF SIMULATED HALOS 

4.i- Overview 

We now study the distributions of simulated halos as a function of Vcirc- The results of 
the computations of A^( Vcirc) A Vcirc scaled to (51.2 Mpc)"^ comoving volumes are presented 
in Figs. 13 through 16. We include observational estimates using both F — 1 and F = 1.1 
to relate cri to V^irc for comparison. 

We focus our efforts on the following three simulations: CDM1(128"^,51.2,280), 
CDM12(64^51.2,52), and CDM16(1443, 100,85). CDMl offers good mass resolution 
(^part = 4.4 X 10^ Mq), CDM12 offers good force resolution (Plummer softening of 
40 kpc comoving), and CDM16 offers fairly good mass and force resolution (mpart = 
2.3 X 10^° Mq] Plummer softening of 65 kpc comoving) yet is computed in a 100 Mpc 
box. Again, we demonstrated in § 2 that the CMF is not very sensitive to the box size 
but we demonstrate in Paper II that clustering statistics require boxes larger than 51.2 
Mpc on a side. 

Our goal in the following sections is to attempt to constrain the amplitude of the 
primeval density fiuctuations of the = 1 CDM model from halo circular velocity distri- 
butions. We devote separate subsections for both high mass and low mass halos, which 
require special treatment for determining reliable simulated and observed distributions. 

4-2. Circular Velocity Distributions of Simulated Halos 
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We measure A'"(Fcirc)^Kirc from the simulations. We ask the questions: 1) Over 
which range of circular velocities do the results agree with the observations? 2) Over 
which range of circular velocities do the results disagree with the observations? 3) Do 
the results depend on numerical resolution and techniques for identifying halos? 

In Fig. 13 we show N (Vdrc) ^Vdrc at as — 0.5, 0.7 and 1.0 for an analysis of 
CDM1(128^51.2,280) using 512^ DENMAX and FOF(Z=0.1) and FOF(Z=0.2). We see 
from Fig. 13 that the number of halos agrees with the observations very well from about 
150kms-i to 350 kms"^ for DENMAX and FOF(;=0.2). The results for FOF(Z=0.1) do 
not fare as well. These statements are true for all three epochs; however, the excess num- 
ber of massive halos gets worse with increasing erg. DENMAX is a compromise between 
FOF(/=0.2) which sometimes merges halos and FOF(/=0.1) which fails to produce some 
halos. These results are encouraging for studies that use FOF(Z=0.2) such as Frenk et 
al. (1988); however, FOF(Z=0.2) occasionally links together visually distinct halos. 

We now study the effects of force resolution, choice of DENMAX grid, and choice 
of R (used to compute Kirc) on A^("l/circ)A"l4irc- We show N{VciTc)^VciTc at erg = 0.5 for 
CDM12(64^,51.2,52) in Fig. 14. The first thing we notice is that the agreement with the 
number of simulated halos with the observations from 150kms~^ to 350kms~^ is even 
better than it is for the low force resolution PM simulation discussed above, particularly 
for Vcirc ~ 200 km s-^ We also see that the results are not very sensitive to the choice of 
R except for the few very massive halos. This is not surprising since most of the circular 
velocity profiles are fiat beyond 150 kpc comoving except for the most massive halos — cf. 
Fig. 11. 

It is encouraging that the results are not very sensitive to the choice of DENMAX 
grid except for the most massive halos and for the 64^ grid. This is not true for the halo 
masses described by the CMF in § 2 — we show later that this is because the different 
DENMAX grids significantly affect peripheral particles beyond the distance R used to 
compute the circular velocities. It is not surprising that the very coarse 64^ grid fails to 
match up to the finer grids. 

We conclude this discussion by testing the sensitivity of the agreement of the number 
of simulated halos with the observations for the different P^M simulations (Fig. 15) and 
then by studying the evolution of the number of the simulated halos (Fig. 16). 

In Fig. 15 we show iV(Fcirc) AKirc for CDM12(643,51.2,52; e = 40 kpc comoving) 
(top panel) and for CDM16(144^, 100,85; e = 65 kpc comoving) (bottom panel), both 
analyzed with a 512^ DENMAX grid. We extract several facts from Fig. 15. First, the 
trend of increasing number of halos with increasing force resolution is verified comparing 
the simulations with e = 40 kpc comoving (top panel) and e = 65 kpc comoving (bottom 
panel), but the differences are small. We found in § 2 that the CMF was higher for higher 
mass resolution simulations and for higher force resolution simulations independently; but 
here force resolution must be dominating because the e = 40 kpc comoving simulation 
has slightly lower mass resolution than the e = 65 kpc comoving simulation, yet still 
produces slightly more halos at a given Kirc- We also found in § 2 that the differences 
in the CMF versus mass resolution were much smaller when we imposed a radius cut on 
the masses. This is equivalent to computing circular velocities. 
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In Fig. 16 we show Kjirc for CDM16(144^,100,85) using a fixed physical radius. We 
hst the epochs as redshifts, z = 1/a — 1, because here we are studying the evolution of the 
halos for a fixed normalization. We assume that the present epoch, 2 = 0, is ao = erg = 1. 
We keep the physical radius cut constant at 100 kpc by using a lOOao/a kpc comoving 
radius cut in equation (3.2). The vertical axis is scaled to a (51.2 Mpc)^ comoving box, 
however. 

We see in the panels of Fig. 16 that N (Vcirc) ^Vdrc (using a fixed proper radius) 
evolves strongly with redshift. At z — 9.9 halos are still forming. The major era when 
galaxies begin to take on the observed distribution is around z = 3.7 to 2.2. Further 
evolution indicates that the halos are merging, i.e. the curves are decreasing. At inter- 
mediate circular velocities (Kirc = 200kms~^) the number of halos decreases by a factor 
of 3.7 from the maximum at 2; ~ 3.7 to 2; = 0. For smaller halos (Fcirc = 150 km s"-*^) the 
effect is higher, a factor of 4.7 from the maximum at 2~3.7to^ = 0. From z = 0.4 to 
z = the factor is rougly constant at ~ 1.5 over wide range of V^irc- The most massive 
halos grow at the expense of the smaller ones. For halos with total bound masses ex- 
ceeding 2.3 X lO^^M© in CDM16(1443, 100,85) {i.e. 1000 particles), we find 245 halos at 
2 = 1, 292 at 2; = 0.4, and 285 at 2 = 0. Therefore, the number of massive halos, unlike 
the lower mass halos, grows little for 2; < 1. 

The mergers implied by Fig. 16 are interesting in themselves, and they arc impor- 
tant for Paper II where merging forms massive systems which have a profound effect 
on galaxy clustering and velocity statistics. Frenk et al. (1988) also found merging in 
their simulations with decreasing redshift. There exists some observational evidence for 
merging. Excess counts of faint galaxies (Tyson 1988) compared with present galaxy 
populations suggest the possibility of merging (Guiderdoni & Rocca-Volmerange 1990; 
Cowie, Songaila, & Hu 1991). The merger hypothesis is not without controversy, how- 
ever, for other possibilities and complications, highlighted by various authors, include 
1) luminosities may evolve more rapidly for faint galaxies than for bright galaxies (e.g. 
Broadhurst, Ellis, & Shanks 1988); 2) the geometry of the universe may be different from 
Einstein-de Sitter {e.g. Fukugita et al. 1990); or 3) the faint galaxies may represent a 
separate population {e.g. Efstathiou et al. 1991). 

Cowie, Songaila, & Hu (1991) argue that the faint galaxy excess is a factor ~ 4 — 5 
from 2 ~ 0.25 to 2 = assuming no luminosity evolution for these modest redshifts. 
Although CDM16 predicts only a factor of ~ 1.5 from 2; = 0.4 to 2; = 0, we cannot 
accurately address galaxy merging with our dark simulations for the following reasons: 
1) we underestimate merging by always associating one galaxy per halo and 2) wc over- 
estimate merging by always assuming that when halos merge their associated galaxies 
merge. Complications aside, since the reduction is ~ 1.5 over a wide range in Kirc: we 
may naively assume that only ~ (1 — 1/1.5) x 0.7 = 23% of the spirals (assuming a 70% 
spiral fraction) have not experienced a major merger since 2 ~ 0.4. This is problematic 
since Toth & Ostriker (1992) argue that high merger rates in the last 5 Gyr (z=0.37 for 
= 1, Hq = 50 kms~^ Mpc~^) can heat disk galaxies beyond observed levels. Further- 
more, if we were to identify gq — as — 0.5 as the present day. Fig. 16 would still apply if 
the Fcirc values were all multiplied by This would move the ug = 0.5 curve {z = 1.0 

in the figure) into agreement with the Schechter function, as it should from Fig. 15. From 
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this we conclude that in the CDM model merging should continue into the future at a 
rate as prodigious as the recent past further violating the Toth & Ostriker (1992) limits. 
For a more detailed examination of merging in CDM models, see Kauffmann & White 
(1993). 

4.3. Massive Halos: Computational Issues 

Since the number of halos from CDM16(144^,100,85) agrees with the observed num- 
ber of halos in the range 150kms~^ ^ V^irc ~ 350 km s"-*^ we now focus on the discrep- 
ancies outside these ranges. In this section we explore circular velocities at various radii 
and we investigate the sensitivity of the formation of massive halos to dynamic range and 
to methods for identifying the halos. In § 4.4 we compare the number of simulated halos 
with high ar to the number of observed bright ellipticals, followed by a discussion of low 
mass halos in § 4.5. The purpose of this section is to reveal which computational effects, 
and why, affect the massive halos. 

We present the four most massive halos at erg = 0.5 from CDM12(64^,51.2,52) in 
Table 2 and from CDM1(1283,51.2,280) in Table 3. The halos are labeled A, B, C, 
and D. These two simulations use equivalent initial conditions. Corresponding halos are 
identified. In the tables we list the circular velocities in kms"-*^ using R = 150 kpc 
comoving, 200 kpc comoving, and 300 kpc comoving. We also list the bound masses 
[R < 00) in solar masses. In Table 2 the results are tabulated for a 512^, 256^, 128^, 
and 64^ DENMAX grid, all at dg = 0.5. In each column we also list a local rank. The 
number n means the halo is the n^^ largest halo in the catalog using the method for halo 
identification mentioned in the first column. Note that the circular velocity profiles for 
these massive halos are still rising far beyond the softening scale. Here we are interested 
in their profiles at large radii. We use Uj. and V^irc extrapolated to more reasonable radii 
in the next section. 

We use Table 2 to study the effect of the choice of R and the DENMAX grid on 
the massive halos. The first important feature brought out is that Fcirc increases with 
increasing radii. These massive halos have extended halos with rising circular velocities 
at these scales (cf. Fig. 11). The next trend we observe is that the circular velocities, 
unlike the CMF without a radius cut, are not very sensitive to the choice of DENMAX 
grid. However, the slight differences are explained below. 

In Fig. 17 the bound particles from halo B found in CDM12(643,51.2,52) are shown 
using the various DENMAX grids. We see that the coarser DENMAX grids (< 256^) 
merge the massive halo with an additional small halo (located at x 200 kpc comoving, 
y —200 kpc comoving). The mass of this "appendage" is small and is far enough 
away from the core (about 300 kpc comoving) so that it contributes little to the circular 
velocity defined within 300 kpc comoving. Nevertheless, it reveals substructure present 
in the higher force resolution simulation. 

The lower resolution DENMAX grids also lead to the inclusion of more peripheral 
(distant) particles. This is not serious since this does not involve a lot of mass and only 
involves particles well beyond 300 kpc comoving from the halo core. When the DENMAX 
grid is finer than the interparticle separation in the periphery, the density gradients are 
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not present to move the particles into the halo. This fact partially explains why the CMF 
(in § 2), based on total bound masses, is more sensitive than Vdrc to variations of the 
DENMAX grid. 

We now consider the effects of force resolution. In Fig. 18 we show the same halo B 
but from the low force resolution PM simulation CDM1(128^,51.2,280) — we show every 
eighth particle for comparison with the 64^ particle P^M simulation. The force resolution 
is too low to produce the "appendage" that we see in the P^M simulation — therefore, 
there is no significant difference between the 512^ grid DENMAX and the 256^ grid 
DENMAX results. We conclude that high force resolution reveals more substructure 
than low force resolution and that high resolution DENMAX grids are required to reveal 
this substructure. 

We see in the lower right panel of Fig. 18 a major failing of FOF(/=0.2). This is 
a particularly pathological example. Of course we could naturally prune this halo into 
separate halos. It is not practical, however, to examine visually and prune manually the 
thousands of halos produced in each simulation. 

The CDMl PM halos A through D, corresponding to the halos studied in the CDM12 
P"^M simulation, are tabulated in Table 3. The results are shown for two DENMAX grids 
and two FOF linking parameters, all at as — 0.5. From Table 3 we conclude: 1) The 
512^ DENMAX results compare well with the 256^ DENMAX results. 2) The FOF 
analyses fail to agree with the DENMAX analyses. The difference between FOF(Z=0.1) 
and FOF(/=0.2) is not too great since these massive PM halos do not have a lot of sub- 
structure. The exception is Halo B in the FOF(/=0.2) analysis. This is the pathological 
halo shown in the lower right panel of Fig. 18. Three visually distinct halos are merged 
together and the center-of-mass is such that the "halo" is highly non-spherical leading to 
unreliable circular velocities. 

We also find, from Table 3, that the P^M halos are more compact than the PM 
halos. If we compare Vdrc defined at 150 kpc comoving in Table 2 for P^M CDM12 
with Vcirc defined at 300 kpc comoving in Table 3 for PM CDMl, we find comparable 
values of Kirc- By the time we go out to 300 kpc comoving in the PM simulation we 
pick up enough particles to give the same circular velocity as the P^M simulation using 
150 kpc comoving. This is because we choose values of R to be the radius where most 
of the circular velocities are flat. These radii are directly related to the force resolution. 
However, things do not always work out this nicely for the massive halos that have rising 
circular velocity profiles, as we can see by comparing halo D in Table 2 for P^M CDM12 
using 150 kpc comoving and halo D in Table 3 for PM CDMl using 300 kpc comoving. 
The differences in circular velocities are significant enough to shift some of the massive 
halos into adjacent 25 kms~^ bins. 

We now summarize some effects arising from the computational techniques that in- 
fiuence the number of massive halos in the N (Vdrc) ^Vdrc histograms. 1) The results are 
sensitive to the choice of R used to compute the circular velocities — this is obvious since 
the circular velocity profiles are not fiat for the massive halos. 2) We have shown that 
higher resolution DENMAX grids reveal more substructure in some of the massive halos 
found in the higher force resolution simulations. However, from the images it appears 
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that no obvious substructure is present in many of the massive halos. 3) Lower resolu- 
tion DENMAX grids include more peripheral particles in the halos than higher resolution 
DENMAX grids. This arbitrary choice of DENMAX grid does not affect most compu- 
tations of circular velocities. It does, however, affect the computations of total bound 
masses; this explains why the CMF is more sensitive to the limitations of the current 
version of DENMAX than is the case for the circular velocities. 4) FOF(/=0.2) occasion- 
ally links together visually distinct halos. FOF(/=0.1) and FOF(/=0.2) produce similar 
results for many of the massive halos but they often fail to match up with DENMAX 
results which, visually, appear to do a good job in many cases. 5) The P^M simulations 
produce halos that are more compact than the PM simulations. However, if a larger 
value of R is chosen for the PM simulations, then the PM circular velocities agree with 
the P^M circular velocities in most cases. 

4-. 4- Massive Halos: Simulations versus Observations 

The distribution of simulated halos with circular velocities in the range 
150 kms"-*^ ^ Vcivc ~ 350 kms"-*^ is in reasonable agreement with observations (Fig. 15). 
However, there are too many halos with circular velocities exceeding 350kms~^. A sim- 
ulation with increased force resolution can reveal more substructure in massive halos and 
a continuum- limit DENMAX algorithm would be helpful for analyzing such simulations. 
We take the approach, in this paper, that these massive dark matter halos represent 
single, large galaxies. The possibility that they may represent clusters is studied in detail 
in Paper H. 

We use fairly complete catalogs of observed bright ellipticals to estimate their num- 
ber density. It is not accurate enough to estimate the brightest, relatively few elliptical 
galaxies simply from a Schechter luminosity function and a Faber-Jackson relationship. 
The problem is exacerbated by the large amount of scatter relating ui to Fcirc for the 
simulated massive halos. In this subsection we instead use ar to characterize the simu- 
lated massive halos and we compare them with the number of observed ellipticals using 
complete elliptical surveys. We use these comparisons to constrain the normalization of 
the = 1 CDM power spectrum using the fact that as the simulations evolve merging 
creates more massive halos. 

We begin by noting that there are observed galaxy halos that have large measured 
circular velocities beyond ~ 100 kpc. The giant elliptical galaxy M87 has been studied by 
many workers using the X-ray emitting gas to trace the gravitational potential well {e.g. 
Fabricant & Gorenstein 1983; White & Sarazin 1988; Tsai 1994). Tsai (1994) modeled 
the X-ray emission from M87 using a multi-phase gas assumed to be in hydrostatic 
equilibrium. Tsai found the best fit gas temperature and mass density profile consistent 
with both X-ray continuum and line emission data. His results are consistent with the 
velocity dispersions of Sargent et al. (1978) and Mould et al. (1990) on small scales. 
(Note that the mean, radial velocity dispersion, a^, of stars in M87 from 1 kpc to 4.5 kpc 
is roughly constant at only 278±11 kms~^, yet can be as high as 350kms~^ well within 1 
kpc-cf. Sargent et al.) The inferred mass within 300 kpc assuming a Hubble constant of 
50 kms"-*^ Mpc~^, is approximately 2.5 x lO^^M© with a corresponding circular velocity 
of 592 kms"-^. Curiously, 1.1 x 592 kms~-^/"\/3 = 376kms~-^ (see equation (3.8)) which is 
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close to the 350 km s"-*^ measurement (within 1 kpc) from Sargent et al (1978). However, 
since it is not clear which small-scale star measurements should be related to large-scale 
dark matter measurements, we adopt an empirical scaling law which relates Faber et al. 
(1989) central velocity measurements (used as a complete catalog of nearby ellipticals) 
to Tsai (1994) circular velocity measurements on large scales. 

In Fig. 19 we show circular velocity profiles for halos B and C (see Tables 2 and 3) 
from the simulations at erg = 0.5 and for M87 (Tsai 1994). We choose halo C because it 
has a circular velocity comparable to M87 at large radii. (Halo B has a higher circular 
velocity than M87.) The profiles from PM CDMl rise slowly which is expected since the 
force softening is 280 kpc comoving. The profiles from P^M CDM12 rise more quickly 
than PM CDMl because of higher force resolution. Ignoring the fact that many of the 
simulated halos are still rising beyond 150 kpc comoving, the conjecture that at least 
some of the very massive simulated halos are similar to objects like M87 is seen to be 
plausible. 

We offer a possible explanation why the simulated rotation curves are still rising be- 
yond the softening scale for P^M CDM12 while M87 has a very fiat rotation curve. During 
the dissipational formation of M87, dark matter can be pulled into the central region by 
baryonic infall {e.g. Blumenthal et al. 1986). If we examine Vdrdi^) in figure 3 from 
Blumenthal et al. (though from a system with maximum circular velocity ~ 200kms~^) 
we estimate that the ratio of the distances where Vdrc turns over is ~ 80 kpc/10 kpc = 8. 
For the P^M CDM12 halo C profile shown in Fig. 19, this effect could possibly "pull" 
the turnover in the dark matter rotation curve from ~ 200 kpc to ~ 25 kpc, consistent 
with the turnover in the M87 profile shown in Fig. 19. 

We examine the largest halos found in the simulations and we compare them to 
one of the most massive and luminous galaxies known — the central cD galaxy in the 
cluster A2029 (Dressier 1979; Uson, Boughn, & Kuhn 1991). The mass profile of this 
galaxy has been estimated with a 3-component model by Dressier (1979): 1) a "normal" 
elliptical galaxy, 2) an extended halo of luminous material out to 100 kpc, and 3) a 
dark cluster-filling component. Dressier estimated the mass within 100 kpc (for Hq = 
50kms-^ Mpc"^) to be ~ 3.9 x 10^^ M© and within 1 Mpc to be ~ 8.3 x 10^^ Mq. 
The evidence that the material within 100 kpc is part of the central cD galaxy is strong, 
but there is some controversy about the mass out to 1 Mpc. Uson, Boughn, & Kuhn 
have argued that the material out to 1 Mpc and beyond is indeed part of the central 
cD galaxy. They measure diffuse light out to several Mpc. They found that it has an 
elliptical profile with the same axis ratio and orientation as the central cD galaxy, and 
that this is different from the distribution of the cluster galaxies as a whole. 

To compare with the above measurements, we compute the mass within 100 kpc 
comoving and 1 Mpc comoving from the simulated halos. Using CDM12(64^,51.2,52) 
we find the halo with the largest mass within 100 kpc comoving and the halo with the 
largest mass within 1 Mpc comoving. The results within 100 kpc comoving (more than 
twice the Plummer softening) are 9.0 x at as = 0.5, 1.7 x lO^^M© at as = 0.7, 

and 2.9 x IO^^Mq at ug = 1.0. None of these are greater than Dressler's estimate for the 
central cD galaxy in A2029, 3.9 x lO^^M©. Within 1 Mpc comoving we find 1.3 x lO^^M© 
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at as = 0.5, 1.9 x IO^^Mq at ag = 0.7, and 1.9 x IO^^Mq at erg = 1.0. Again, these are 
all smaller than Dressler's estimate, 8.3 x IO^^Mq. Thus, we cannot rule out CDM by 
arguing that it produces halos with absolutely too much mass. We also cannot rule out 
CDM merely by the fact that our simulation fails to make at least one halo as massive 
as the central cD galaxy in A2029 — we sample only a 51.2 Mpc box while A2029 is at a 
distance of 470 Mpc. 

We examine other simulations in an attempt to find halos as massive as the central cD 
galaxy in A2029. We examine CDM6(256^,51.2,190) at R= 1 Mpc comoving at = 1.0. 
The most massive halo at this radius has a mass of 3.7 x lO^^M©. Also, we examine 
CDM16(144'^, 100,85) at erg = 1.0 with no cut in radius, and the most massive halo has a 
mass of 8.9 x 10^^. This is a larger box with larger waves in the initial conditions and a 
different set of initial random numbers. The model still fails to produce a halo as massive 
as ~ 8.3 X lO^^M©. Thus far the 1] = 1 CDM model may be safe. 

Although we cannot reject CDM based on the most massive halo in the simulations, 
maybe we can reject it based on the large number of slightly less massive halos that are 
formed. Because the most massive galaxies are ellipticals, we compare the number of 
simulated halos with large radial velocity dispersion ar (second of equations (3.1)) with 
the number of ellipticals having large line-of-sight central velocity dispersion. For the 
observations we use the samples of nearby bright elliptical galaxies from Faber et al. 
(1989) and from the Dressier (1991) supergalactic plane redshift survey. We count the 
number of ellipticals in these samples with log^^g ^ 2.5 (20 ellipticals). We then impose 
a distance cut of 6000 km s~^ (based on corrected distances from column 12 of table 3 
from Faber et al. 1989). This leaves 14 ellipticals. The samples are fairly complete. 
For the range in apparent magnitudes of our list of ellipticals, the completeness fraction 
ranges from 100% for Bt ,^11.6 down to 20% for ellipticals in the southern sample with 
Bt ~ 13. If we fold in the completeness fractions (figure 2 from Faber et al. 1989) the 
number of 14 ellipticals with log^g Cr > 2.5 within a distance of 6000 km s~^ might be as 
high as ~ 23. 

We also estimate the number of observed ellipticals from the above samples with > 
350 km s"-*^. This allows us to study very high values of (7^, for which the completeness 
fractions are much higher. Within 6000 km s~^ there are only three known galaxies with 
ar > 350kms-i: SPS 1120 {ar = 382 km s'^; Bt = 12.68), NGC 507 {ar = 366 km s'^; 
Bt = 11.63) and NGC 4486 (M87; = 361 kms'^; Bt = 9.52). The completeness 
fraction (based on Bt) for SPS 1120 is ~ 30% and the completeness fraction for the 
other two objects is 100%. This tightly constrains the number of observed ellipticals 
with ar > 350kms~^ within 6000 km s"-*^ to 5. 

To demonstrate the inaccuracies at the high mass end associated with methods 
presented in § 3, we compare the above complete estimates with the use of the Faber- 
Jackson relationship (equation (3.7)) using ai (first of equations (3.1)) and the Schechter 
function described in § 3. Again, we weight the Schechter function by 30%; i.e. we only 
estimate the elliptical contributions. We find in a spherical volume of radius 120 Mpc 
comoving: 39 objects with ai > 316kms~^ and 11 objects with ai > 350 km s"-*^. These 
numbers are about a factor of 2 larger than the estimates given above, suggesting that 
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our assumed Faber- Jackson relation underestimates the luminosity and/or that ellipticals 
make up less than 30% of bright galaxies. (We combine ellipticals and lenticulars in 
getting the population fraction of 30%, but lenticulars are undcrrepresented among the 
most massive galaxies.) In any case, these results suggests that we have overestimated the 
observed A''(ycirc) in Figs. 13 through 16 for V^irc ^ 500 km s"-*^, making the disagreement 
with the simulations even worse. 

To make a better comparison of the simulations with observations, we estimate 
the number of simulated halos from four simulations with ctj. > 316 km s""*^ and ctj. > 
350kms~^. The results are shown in Table 4 for simulations with a variety of force 
resolution and mass resolution. All numbers are scaled to a (51.2 Mpc)^ comoving volume. 
The observations are shown as OBS. I (without completeness fractions folded in) and as 
OBS. II (with completeness fractions folded in). All simulations use a 51.2 Mpc box 
except for CDM16 which uses a 100 Mpc box. The initial conditions for CDM6, CDMl, 
and CDM12 are all generated from the same set of 256^ random numbers. The initial 
conditions for CDM2-5 and CDM16 are all generated from different sets of random 
numbers. We also show averages, with la fluctuations, computed from CDMl-5. 

We estimate the velocity dispersions from the simulated halos in two ways: 0"^ and 
a which we describe below. We compute the radial velocity dispersion, second of equa- 
tions (3.1), within a radius listed in the footnotes of Table 4. We count the number of 
halos with (t^ exceeding 316kms~^ and 350kms~^. We also try using ai (not shown), 
first of equations (3.1), and the results are similar to the results using 0"^ (the differences 
arise from the fact that cri is typically ~ 20% lower than cr,. as mentioned earlier). 

The high velocity dispersions of the dark matter may not correspond to the ve- 
locity dispersions of the optical galaxies expected to be embedded well within the 
cut-off radii used here. As a crude estimate of the central velocity dispersion of a 
galaxy expected to be embedded in the simulated halo, we use the following: a = 
361kms~^ X ycirc(-R)/600kms~^. We choose this because the central velocity disper- 
sion of M87 is 361kms~^ from Faber et al. (1989) and the estimated circular velocity 
profile for M87 (from Tsai 1994; see Fig. 19) is about 600 kms"^ for R > 100 kpc. Ad- 
mittedly this is a crude estimate, but it is a simple attempt to use a single, well-measured 
object to scale the simulated data, and it serves as a conservative check for our compar- 
isons with observations. (Note, in this case equation (3.8) works fairly well for F = 1.1. 
However, it may not work well for all objects.) 

In Table 4 we compare the results from the simulations with the observations. The 
most important conclusion is that all cases for ug > 0.4 yield far more halos than the 
observed numbers. The epoch erg = 0.3 is not ruled out since it is difficult to make 
conclusions based on zero or one halo. The results at erg = 0.4 indicate that there are 
too many halos with cr^ > 316kms~^. The problem is less severe (yet not an order of 
magnitude different from cr^.) using a, but the observations tell us that even a single halo 
with ar > 350kms~^ in a 51.2 Mpc box is too high by at least a factor of 10. 

White et al. (1987), at ug = 0.4 using the same normalization of the CDM power 
spectrum as we do, found a single halo with a circular velocity exceeding 800 kms~^ 
in a 50 Mpc box from three simulations, corresponding to 0.36 halos for a single 51.2 
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Mpc simulation. Our CDM12 simulation has nearly identical force and mass resolution. 
Our largest halo in CDM12 at erg = 0.4 has a circular velocity of 567 kms~^ defined at 
100 kpc comoving. We cannot safely rule out erg = 0.4 particularly since results using 
a only reveal one very massive halo. The problem becomes rapidly worse for larger cts; 
erg ^ 0.7 predicts more than 20 times too many galaxies with ar > 350kms~^. This is a 
severe problem for $7 = 1 CDM since estimates of cxg based on clustering typically require 
(7g > 0.4 (Davis et al. 1985; Park 1990; Couchman & Carlberg 1992). 

We now compare the different simulations with each other. We see the general trend, 
in Table 4, that both an increase in mass resolution and an increase in force resolution 
increase the production of massive halos. For the simulations using equivalent initial 
conditions (CDM6, CDM12, and CDMl) we examine corresponding massive halos. The 
higher force resolution simulations produce more compact halos than the lower force 
resolution simulations; the cut-off radius is chosen to compensate for this fact for reasons 
discussed earlier. In some cases, however, the cut-off radius does not compensate for the 
compactness of the high force resolution halos. We also find that the high force resolution 
halos have higher central velocity dispersions. 

An increase in mass resolution also increases the production of massive halos. The 
effect is strongest at cxg = 1.0 where the numbers from the 256^ particle simulation are 
far higher than the other simulations compared with earlier epochs. We observe that the 
results from the 100 Mpc box simulation CDM16 with good mass and force resolution are 
in reasonable agreement with the other simulations except at erg = 1. We also observe 
that the 128^ particle, R1/2 = 280 kpc comoving PM simulations produce the smallest 
number of massive halos — these simulations rank low in the combination of force and 
mass resolution and R = 300 kpc comoving is too close to R1/2 — 280 kpc comoving. 

We conclude that the = 1 CDM model is in serious trouble. The simulations 
produce far too many massive halos and an increase in force and mass resolution only 
make matters worse. We are able to rule out all normalizations of the primeval density 
fluctuations with erg ^ 0.5. Using complete catalogs of nearby bright ellipticals, we have 
constrained the CDM model more convincingly than by using the luminosity function at 
the bright end (cf. Fig. 15 and Frenk et al. 1988). The case against erg = 0.4 is not 
as strong as the case against erg ^ 0.5. We found, at the very least, a single halo with 
an estimated central velocity dispersion exceeding 350kms~^ in a single 51.2 Mpc box 
simulation. The observations predict that we should only find one such object in no fewer 
than 11 simulations. 

We know that the simulations suffer from the overmerging of massive halos. Gas 
dynamical dissipation could reduce the merging of galaxies. The result might be to 
prevent the formation of excessively massive galaxies, although we consider this unlikely 
because dissipation should only increase the central concentration of mass in the most 
massive halos, thereby increasing further the central velocity dispersions. Also, if the 
most massive halos actually should represent clusters of galaxies, then these clusters 
must still have the correct multiplicity function (distribution of richness). Bahcall & Cen 
(1992) concluded that the CDM model with ag = 1.05 produces an order of magnitude 
too many rich clusters. In Paper II we investigate the cluster multiplicity function in 
detail using our own high-resolution N-body simulations. 
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4-5. Low Mass Halos 



We now examine the low mass halos. We found earlier that the P^M simulations 
produce too many halos with Vdrc ~ 150 kms"-*^ (see Fig. 15). Frenk et al. (1988) ar- 
gue that the number of halos is in reasonable agreement with the observations down to 
about 60kms~^ using 32000 particle P"^M simulations in 14 Mpc boxes. However, they 
warned the reader that simulations in larger volumes predict too many halos (White 
et al. 1987). The particle mass in the Frenk et al. simulations is 5.8 x IO^Mq and 
the force resolution is e = 14 kpc. We have two P^M simulations with lower mass and 
force resolution that are computed in a 51.2 Mpc box and a 100 Mpc box giving us 
better statistics: CDM12(643,51.2,52; mpart = 3.5 x lO^^M©; e = 40 kpc comoving) and 
CDM16(1443, 100,85; mpart = 2.3 x 10^°Mq; e = 65 kpc comoving). Using these sim- 
ulations, we explore the effects of resolution and we re-examine the observational data 
at the low mass end in order to explore the apparent excess number of low mass halos 
compared with the observations. 

The smallest galaxies for which there are reliable mass estimates have Vdyc down to 
about 50 kms~^ (see Kormendy 1991 and references therein). Halos from the 40 kpc 
comoving Plummer simulation CDM12 with a cut-off radius of 100 kpc comoving (roughly 
twice the Plummer softening) and 5 particles have a circular velocity of 87kms~^. Halos 
from the 65 kpc comoving Plummer simulation CDM16 with a cut-off radius of 150 kpc 
comoving and 5 particles have a circular velocity of TOkms""*^. Therefore, we can only 
study halos down to 70kms~^ using the P^M simulations. 

Halos from the PM simulation CDM1(1282,51.2,280) with 5 particles and a cut-off 
radius of 300 kpc comoving have a circular velocity of 18kms~^, and halos from the 
PM simulation CDM6(2563, 51.2, 190) with 25 particles and a cut-off radius of 300 kpc 
comoving have a circular velocity of 14kms~^. However, these PM simulations have 
poor force resolution. We show that higher force resolution increases the number of low 
mass halos. Therefore it is misleading to compare the number of low mass halos with the 
observations using the PM simulations. 

Another problem stems from the fact that we need to use large cut-off radii to 
characterize the circular velocities in the PM simulations. In order for a galaxy to undergo 
"complete collapse" in a spherical, Q = 1 model, it has to have an over-density exceeding 
5g = 5p/p ~ 170 (Gunn & Gott 1972; BBKS). A simple calculation shows that this places 
a lower limit on the circular velocity for a given cut-off radius R. The circular velocity 
within R for a density p is simply 



If we demand that the over-density exceed 5g = p/pcrit — 1, where pcrit is the density for 
an ri = 1 universe given by 3HQ/(8nG), we get the minimum allowed circular velocity: 




(4.1) 




(4.2) 
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For 6g = 200 (close to the critical value, chosen to yield a simple formula), Hq = 
50kms~^ Mpc~^, and a comoving cut-off radius R measured in kpc, we arrive at the 
simple formula for the minimum allowed circular velocity in kms~^: 



Eq. (4.3) puts a severe limit, Vdrc ~ 150 km s"-*^, on the PM simulations that require 
R = 300 kpc comoving. For the P^M simulations that require R = 100 kpc comoving 
and 150 kpc comoving, the restrictions are 50kms~^ and 75kms~^ respectively. 

Before exploring the simulations, we need to examine the observational parameters 
used for the TuUy-Fisher relationship and the Schechter luminosity function (see § 3) for 
faint galaxies. Since we have already shown that the simulations appear to produce too 
many halos at the low mass end, we conservatively choose parameters that produce the 
largest number of low mass halos allowed within the uncertainties of the observations. (We 
find that there are still too many halos predicted by the CDM model so we are not forcing 
the observations to agree with the model — we are simply estimating how significant is 
the discrepancy.) In the following discussion we re-scale all relevant numbers to a Hubble 
constant Hq = 50kms~-'^ Mpc~^. 

First we consider the TuUy-Fisher relationship in equation (3.6). Pierce & TuUy 
(1988) reported that the scatter in this relationship is ±0.25 magnitudes. In their fits 
(figure 9 from their paper), they found that the faintest galaxy studied, Mbj, ~ —16 -|- 
5 log]^o(50/85) ~ —17.2, is slightly brighter than predicted by their best fit. Alternatively, 
if one measures the circular velocity of this faint galaxy, the Pierce & TuUy relationship 
would predict that the galaxy is fainter than it actually is. Since the luminosity function 
is an increasing function of decreasing luminosity, one would overestimate the number of 
faint galaxies. We take an extreme point of view. We will use equation (3.6) as is with 
an added value of 0.5 magnitudes — this is twice the reported scatter quoted by Pierce & 
TuUy and it results in an increase in the estimate for the number of observed halos as a 
function of Vcirc- 

Next, we consider the luminosity function. The estimates of Efstathiou, Ellis, & 
Peterson (1988) for the parameters of the luminosity function are estimated to hold 
down to about M^j, ~ — 16-1-5 log]^Q(50/100) ~ —17.5. The luminosity function has been 
studied by previous workers down to comparably faint magnitudes (see Felten 1977 for a 
review). This faint limit is comparable to the faint limit of the TuUy-Fisher relationship. 
Therefore, we use the parameters of the Schechter luminosity function given in § 3 but 
we use the reported errors to yield the maximum number of faint galaxies. We assume 
that 100% of the faint galaxies are spirals. We use $* = (1.56 + 0.34) x lO'^/i^ Mpc"^ 
and M^^ = (-19.68 - 0.10) - 2.5log^oh-^ with h = 1/2 and we use a = -1.07 - 0.05. 

These changes in the TuUy-Fisher relationship and the Schechter luminosity function 
increase the estimated number of faint halos with 50kms~^ < Vcirc < 75kms~^ from 
373 to 582 galaxies in a (51.2 Mpc)"^ comoving volume. When we show the number of 
observed low mass halos, we use both the parameters given in § 3 and the extremely 
stretched parameters given in this section. 




kms 



-1 



(4.3) 
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In Figs. 20 and 21 we show the low mass end, 50kms~^ < Vdrc < 200kms~^, from 
the simulations. We still use 25kms~^ wide bins but we re-bin the data from 50 kms~^ 
to 75 kms~^, 75 kms~^ to 100 kms~^, etc. The observations using the parameters 
described in § 3 with F — 1 are shown as solid squares. The "maximum" number of faint 
halos allowed by the observations minus the default values is used for the ± error bars 
(the asymmetry is because we use logarithms on the vertical axes; note these are not la- 
error bars). 

The figures list the various simulation parameters and the choices of R. The simu- 
lations are shown down to circular velocities such that the bins are complete given the 
mass resolution limit. These restrictions exceed the restrictions based on the over-density 
argument, equation (4.3), for the P^M simulations. The PM simulations are restricted by 
the over-density argument to VcWc ~ 150kms~-^. The PM simulations produce fewer low 
mass halos than the P^M simulations. This must not be taken to mean better agreement; 
instead it is an example of how poor force resolution can give misleading results. 

The results for the e = 40 kpc comoving Plummer simulation, CDM12, and the 
e = 65 kpc comoving Plummer simulation, CDM16, are in reasonable agreement with 
each other above 100 kms~^. CDM16 has slightly more power on small scales, ANyquist = 
W^Nyquist = 2 X (100 Mpc/144), than does CDM12, ANyquist = 2 x (51.2 Mpc/64). We 
learned from the CMF studies (§ 2.3) that small-scale waves affect the low mass end. 

In Fig. 21 we show results from CDM12 for R = 100 kpc comoving and 150 kpc 
comoving at ag ~ 0.3, 0.4, 0.5, 0.7, and 1.0. We see that the results are not very sensitive 
to R. In all cases, there are still too many halos particularly below 125 kms~^ and 
definitely below lOOkms"^. We see that the number of low mass halos, unlike the high 
mass halos, decreases with increasing expansion factor (both effects are due to merging) . 

We now compare a few numbers at as — 0.4 and erg = 1.0 from CDM12 (with 
R = 150 kpc comoving) for the ranges 75kms~-^ < VcWc < lOOkms"^ and lOOkms"-^ < 
V^irc < 125kms~^. The numbers in these bins from the simulation are 1087 and 495 
respectively for erg = 0.4 and 724 and 333 respectively for erg = 1.0. Using the ob- 
servational parameters from § 3 we find 240 and 168 respectively. Using the extreme 
observational parameters discussed in this section boosts the numbers to 360 and 247 re- 
spectively. Therefore, the excess number of halos below VcWc ~ 125kms~-^ is significant. 
The simulations produce factors ~ 2 — 3 too many faint halos. As a final check, we use 
the parameters described in this section but we try a = —1.25 which boosts the number 
for 75kms-^ < Kirc < 100 km s"^ to 553— stiU short of the 724 to 1087 found in the 
simulation. 

We conclude that the Q — 1 CDM model produces too many low mass halos com- 
pared with the observations for VcWc ~ 125kms~^. We have compared the numbers from 
a 40 kpc comoving Plummer simulation with the largest estimates allowed by the obser- 
vations and the discrepancy is still large (about a factor of 2). Increased force resolution 
and increased small-scale power in the initial conditions make the disagreement worse. 
Although these disparities are large, Dekel & Silk (1986) argued that supernovae in dwarf 
galaxies can cause significant gas loss, and therefore dim the galaxies with small Fcirc- 
Perhaps the TuUy-Fisher relation breaks down at such small Vdrc (though there is little 
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indication of this in the data of Pierce & TuUy 1988). For these reasons, though, we 
consider the excessive number of low mass halos in the CDM model to be less serious 
than the excessive number of high mass halos. 

5. CONCLUSIONS 

A promising result for the CDM model is that the distribution of halos as a function 
of circular velocity agrees rather well with the observations for circular velocities in the 
range 150kms~^ to 350 km s^^. The agreement is better over this range for the P'^M 
simulations versus the lower force resolution PM simulations and the agreement is not 
very sensitive to a Plummer softening of 40 kpc comoving versus 65 kpc comoving over 
this range. However, we found serious problems outside of this range and the problems 
are made worse by increasing the force resolution and the mass resolution. Although 
CDM16 is not the highest resolution simulation, it is computed in a 100 Mpc box; we will 
discover in Paper II that 51.2 Mpc boxes are too small to accurately study clustering. On 
the other hand, the properties of individual halos are not very sensitive to the differences 
between a 51.2 Mpc box and a 100 Mpc box — this fortunate fact allowed us to use many 
of the 51.2 Mpc box simulations to explore effects arising from varying mass and force 
resolution and from different methods for identifying halos. 

We now summarize the chief conclusions found in the preceding sections. 

1. We studied the cumulative mass fraction CMF(M), the fraction of all the mass in 
halos more massive than M. We found the following: (a) We need to compare the CMF 
from simulations analyzed with the same effective DENMAX resolution — lower resolution 
grids include more peripheral particles, increasing the total masses of the halos; (b) The 
simulation-to-simulation scatter is small except for the most massive halos; (c) Higher 
mass resolution and higher force resolution each increase the CMF independently. The 
effect of increased mass resolution on the CMF is reduced if we impose a distance cut 
from the density peaks of the halos, and comparisons of simulations with 64^, 128^, and 
256^ particles indicate that the convergence of the CMF with such a cut is plausible, (d) 
Small-scale waves in the initial conditions have a very small effect on the CMF except for 
the smallest halos; (e) Long waves (with wavelength exceeding 51.2 Mpc comoving) in 
the initial conditions do not affect the CMF for amplitude a% < 0.5; and (f) The Press- 
Schechter theory with 5c = 1-68 predicts too many massive halos and a more rapid growth 
of the CMF than found in the simulations. Substructure within halos is apparently not 
erased as rapidly as implied in the Press-Schechter theory. 

2. Simulated halos generally have mass distributions characterized by flat rotation 
curves extending from about two softening radii to 500 kpc comoving or more. The 
most massive halos have shallower density profiles, resulting in rising rotation curves. 
The independence of circular velocity with radius for most halos allows us to compare 
simulated halos at radii of 150-200 kpc comoving (in the P^M simulations) with real 
spirals at 10 kpc comoving or less. 

3. The distribution of circular velocities of simulated halos was compared with 
observations. We noted above the good agreement for 150kms~-^ ^ VcWc ~ 350kms~^ 
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for any of the three normahzations as = 0.5, 0.7, and 1.0. In the analysis, for this range in 
circular velocities, we found the foUowing: (a) The agreement with the observations is best 
for the P"^M simulations and is not very sensitive to simulations with a Plummer softening 
of 40 kpc comoving versus 65 kpc comoving. (b) The results from DENMAX agree better 
with FOF(Z=0.2) than with FOF(Z=0.1). (c) The distribution of circular velocities, unlike 
the CMF, is not very sensitive to the DENMAX grid, but higher resolution grids are 
required to pick out substructure in the P^M simulations, (d) The number of halos 
characterized by their circular velocities (using a fixed, physical radius) indicates, if erg = 1 
is the present epoch, that the galaxy mass function takes on its present shape by ^ ~ 3.7. 
Between this epoch and z = 0, merging reduces the number of halos by about a factor of 
3.7. Merging is predicted to continue into the future. 

4. We conclude from the studies of massive halos that the Q = 1 CDM model is 
in trouble if these systems represent individual galaxy halos. We are able to rule out 
normalizations of the primeval density fluctuations with erg ^ 0.4 based on the number of 
massive halos if the halos represent individual galaxies, although the lower limit for cxg is 
uncertain. We compared the simulations not only with the observed luminosity function, 
but also with complete samples of bright nearby ellipticals. These observations constrain 
the model to erg < 0.5. We cannot rule out CDM based on the most massive halo — we do 
not find any halos at any epochs with masses exceeding the inferred mass of the central 
cD galaxy in A2029. If the massive halos represent unresolved clusters, with the central 
galaxy having a smaller central velocity dispersion than the surrounding halo, we may 
relax these constraints. We consider this possibility further in Paper II. 

5. We conclude from the low mass studies that the f2 = 1 CDM model pro- 
duces too many low mass halos (by factors ~ 2 — 3) compared with the observations 
for V^irc ^ 125kms~^. The number of faint halos decreases with increasing erg because 
of merging. Nevertheless, the excess is significant even at erg = 1 using extreme assump- 
tions about the observational uncertainties. We do not find reasonable agreement down 
to ~ GOkms"-*^ as reported by Frenk et al. (1988). Gas loss in dwarf galaxies (Dekel 
Sz Silk 1986), however, might dim a significant number of dwarf galaxies, making this 
problem less critical for CDM than the high mass problem. 
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FIGURE CAPTIONS: 



FIG. 1: Cumulative mass fractions for CDM1(128^,51.2,280) analyzed using 
FOF(/=0.1), FOF(/=0.2), and 512^ DENMAX. a) Compares DENMAX (solid curves) 
with FOF (dot-dashed curves for /=0.1; dashed curves for Z=0.2) and b) the effect of the 
removal of unbound particles (solid curves for bound particles; dot-dashed curves for all 
particles). Each case has three curves — lower curves (cxg = 0.5), middle curves (cxg = 0.7), 
and upper curves (erg = 1). 

FIG. 2: Cumulative mass fractions for 512^ DENMAX halos from CDMl-5. All use 
128^ particles, a 51.2 Mpc box, and a force softening distance R\/2 = 280 kpc comoving. 

FIG. 3: Cumulative mass fractions for 512^ DENMAX halos from various simulations in 
51.2 Mpc boxes. The effects of particle number, and force softening, -R1/2, are shown. 
All four simulations are generated from an equivalent set of random numbers. 

FIG. 4: Cumulative mass fractions at erg = 0.5 for DENMAX halos from 

CDM12(64^,51.2,52) using a 512^ grid (solid curve), and lower resolution DENMAX 
grids: 256"^ grid (short-dashed curve); 128"^ grid (long-dashed curve); and 64^ grid (dot- 
dashed curve). 

FIG. 5: Cumulative mass fractions for 512^ DENMAX halos from two simulations in big- 
ger boxes: CDM11(1283, 102.4,560; solid curves) and CDM16(1443, 100,85; short-dashed 
curves). Also, top panel (ag = 0.5) for 256^ DENMAX halos from CDM12(643, 51.2,52; 
long-dashed curve). 

FIG. 6: Cumulative mass fractions for 512^ DENMAX halos from three R1/2 = 280 kpc 
comoving PM simulations in 51.2 Mpc boxes. All three simulations use equivalent initial 
conditions. The dashed curves are for 128^ particles — hut the initial displacements were 
interpolated from the 64^ particle case (dot-dashed curves). 

FIG. 7: Cumulative mass fractions for 512^ DENMAX halos at (Tg = 0.5 from three 
PM simulations using equivalent initial conditions. We include only particles within 
300 kpc comoving from the density peak, a) Compares three simulations indicated by 
particle number N and force softening -Ri/2. The halos whose raw mass (no cut in 
radius and no unbinding) exceeding the transition mass (1.1 x lO-'^^M©; vertical bar) 
have not had their unbound particles removed, b) Results with and without the removal 
of unbound particles above the transition mass, c) The same curves from Fig. 7a except 
the 256^ particle, i?i/2 = 190 kpc comoving simulation has been scaled to R1/2 — 280 
kpc comoving. 

FIG. 8: Cumulative mass fractions from CDM12(643,51.2,52) and CDM16(1443, 100,85) 
and the Press-Schechter theory (PS), all for ug = 0.5, 0.7, and 1 with the CMF being 
larger with increasing erg. a) Dotted curves are PS with a top hat window function 
and 5c = 2. Solid circles are CMFs of raw masses for CDM16 and FOF(/ = 0.2) while 
crosses are CMFs of raw masses for CDM12 and FOF(/ = 0.2). b) Same as a) except 
PS is for a gaussian window function with 5c — '2. and the simulations are CMFs of 
raw masses analyzed with FOF(Z = 0.1). c) PS is for a gaussian window function with 
5c = 2. Solid circles are CMFs from CDM12 computed with raw masses using a 512^ 
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grid DENMAX analysis. Crosses are CMFs also from CDM12 and also computed using 
a 512^ grid DENMAX analysis — however, only bound particles within 200 kpc comoving 
of the DENMAX peak are used to compute the CMFs. 

FIG. 9: cTi and V^rc = {GM/Rf/'^ from CDM16(1443, 100,85) at as = 0.7. Each point 
represents one halo. We show various comoving cuts, a) (Ti(200 kpc) versus cri(100 kpc). 
b) ^^€1^(200 kpc) versus K:irc(100 kpc). c) Fcirc(300 kpc) versus ycirc(200 kpc). 

FIG. 10: V"circ/c"i versus Vchc (all computed using R = 200 kpc comoving) from 
CDM16(1443, 100,85). SoHd Hues are for Vcirc/f^i = ^3, or F = 1 in equation (3.8), 
and dashed lines are for Vdrc/o"! = \/3/l.l, or F = 1.1 in equation (3.8). 

FIG. 11: Circular velocity profiles for 512^ DENMAX halos from CDM12(643,51.2,52; 
e = 40 kpc comoving). For each halo we computed Fcirc at 150 kpc comoving: V150. We 
then sorted the halos from large to small l^so- We show the top ten halos and then every 
twentieth halo thereafter, all the way down to V150 = 150 km s""*^. (This procedure is 
done independently at each epoch.) 

FIG. 12: Cumulative velocity dispersion profiles ai{R) for the same halos shown in 
Fig. 11. 

FIG. 13: Distribution function of circular velocity for CDM1(128^,51.2,280) with the 
circular velocities measured at 300 kpc comoving. The results are scaled to a comoving 
volume of (51.2 Mpc)^ in all of the distribution plots (Figs. 13 through 16) for comparison. 
The dot-dashed curves in these plots are for a Schechter function (F = 1 and F = 1.1; 
the latter gives slightly fewer numbers for bright elliptical halos). We use 512^ DENMAX 
(solid histograms), FOF (/=0.1, short-dashed histograms), and FOF (Z=0.2, long-dashed 
histograms). The histograms, high to low values at Vdrc ~ 200kms~^, are for DENMAX, 
then FOF {1—0.2), then FOF (/=0.1). For ag — 1, each method found two halos for 
Kirc > SOOkms"-^. The bins at Vdrc = 875kms~^ and 1050 km s~-^ each contain one FOF 
{1=0.1) halo and one FOF {1=0.2) halo. The bins at Kirc = 925 km s-^ and 1075 km s-^ 
each contain one DENMAX halo. 

FIG. 14: Distribution function of circular velocity for CDM12(64^,51.2,52) analyzed at 
cTg = 0.5 with DENMAX grids of 512^ (solid histograms), 256^ (short-dashed histograms), 
128^ (long-dashed histograms), and 64"^ (dot-dashed histograms). The comoving radii 
used to define the circular velocities are: a) 150 kpc, b) 200 kpc, and c) 300 kpc. The 
512^, 256^, and 128^ DENMAX grids are nearly indistinguishable except for large Vdrc- 
The coarse 64^ DENMAX grid fails to match up to the other histograms. In the bottom 
panel, each grid identified two halos above 700kms~^. The 800kms~^ bin contains a 64^ 
grid halo and a 512*^ grid halo. The 825kms~^ bin contains a 128^ grid halo and a 256^ 
grid halo. The 875kms~^ bin contains a 512"^ grid halo. The 900kms~^ bin contains a 
64^ grid halo, a 128"^ grid halo, and a 256^ grid halo. 

FIG. 15: Distribution function of circular velocity for 512"^ DENMAX halos from a) 
CDM12(643,51.2,52) and from b) CDM16(1443, 100,85). We use a 150 kpc comoving 
distance to compute Vdrc for CDM12 and a 200 kpc comoving distance to compute Vdrc 
for CDM16. The results are shown at os = 0.5 (solid histograms), 0.7 (short-dashed 
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histograms), and 1.0 (long-dashed histograms); they overlap except for high V^irc where 
there are more halos for larger values of erg. 

FIG. 16: Same as Fig. 15b, except that we use a fixed physical radius of 100 kpc to 

compute circular velocities of halos in a (51.2Mpc)^ comoving volume. The epochs are 
listed in terms of the redshift 2; = 1/a — 1 where we take ao = (Js = 1 to be the present 
day. This plot depicts evolution. 

FIG. 17: The bound particles of Halo B (see Table 2) from CDM12(643,51.2,52) at 
as — 0.5 are shown as found by the various DENMAX analyses. The images are shown 
as X — y projections in units of comoving kpc. The 512^ DENMAX analysis is able to 
resolve the small halo present in the other panels (located at x ^ 200 kpc and y —200 
kpc). 

FIG. 18: The bound particles of Halo B (see Table 3) from CDM1(128^, 51.2,280) at 
as = 0.5 found by the 256^ and 512^ DENMAX analyses and by the FOF(/=0.1) and 
FOF(/=0.2) analyses. (For EOF we do not remove the unbound particles.) We show 
every eighth particle to facilitate a comparison with the 64^ particle simulation shown 
in Fig. 17. There is not much difference in the two DENMAX analyses (apart from the 
peripheral particles) because the PM forces are computed on a 256^ grid. The FOF(/=0.2) 
analysis reveals a dramatic shortcoming of EOF — namely the linking together of several 
dynamically distinct halos. 

FIG. 19: The circular velocity profiles at erg = 0.5 for halos B and C (halo B has 
a larger V^rc than halo C) from CDMl (solid curves: 128^,51.2,280; see Table 3) and 
from CDM12(dotted curves: 64^,51.2,52; see Table 2). We show the profile for M87 as 
computed by Tsai (1994) based on analysis of X-ray emission (short-dashed curve). 

FIG. 20: Distributions of simulated low mass halos from four simulations with vari- 
ous comoving cuts R: CDM1(128^,51.2,280, R = 300 kpc; dot-long-dashed histograms), 
CDM6(256^51.2,190, R = 300 kpc; short-dashed histograms), CDM16(1443, 100,85, 
R = 150 kpc; long-dashed histograms), and CDM12(64^51.2,52, R = 100 kpc; dot-short- 
dashed histograms). We estimate the observed numbers using parameters presented in § 
3 (solid squares) with plus/minus error bars (these are "extreme" systematic errors, not 
la error bars). 

FIG. 21: Distributions of simulated low mass halos from CDM12(64^51.2,52). (The 
observed numbers are solid squares with error bars, see Fig. 20.) The results are shown for 
a) R =100 kpc comoving and b) R —150 kpc comoving. The results are shown at erg = 0.3 
(dotted histograms), 0.4 (short-dashed histograms), 0.5 (long-dashed histograms), 0.7 
(dot-short-dashed histograms), and 1.0 (dot- long-dashed histograms). 
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